\section{Introduction}
\label{section:Introduction}
Fourier Neural Operator \cite{li2020fourier} and Deep Operator Network \cite{lu2019deeponet} are two of the most well-developed neural operators to date. Benchmarks indicate that PDE solvers based on them can be superior to classical techniques at least in the regimes of low accuracy \cite{wang2021learning}, \cite{kovachki2021neural}, \cite{grady2022towards}, \cite{pathak2022fourcastnet}, \cite{yang2022generic}. Architectures of FNO and DeepONet are quite different, but they both use sampling to approximate input function and neural networks to approximate the output. We have reasons to believe that these choices are unfortunate when PDEs are considered.

First, the output of the neural operator is a neural network, i.e., essentially a black-box function. In classical numerical methods such as FEM \cite{ciarlet2002finite}, spectral methods \cite{boyd2001chebyshev} and others, the parametrization allows for extracting bounds on function,  its derivatives, and any other local or global information in a controlled transparent manner. For neural networks, this is not the case. Consider neural networks with $\text{ReLU}$ activation functions in $1D$. The network is a piecewise linear function, so on a given interval it is uniquely fixed by partition and the value of derivative on each subinterval \cite{devore2021neural}. The number of subintervals is exponential in the number of layers, so the partition is not easy to find. Clearly, simple strategies such as sampling are not going to work. The example of this situation is given on \cref{fig:bad_network}. Here on grids with $h$, $h/2$, $h/4$, $\dots$, neural network approximates $2x$, and on the grid with $h/2^{n}$, it significantly diverges from $2x$. Such behavior is clearly undesirable from the perspective of scientific computing.\footnote{Perhaps the only argument for neural networks is to use them for high-dimensional problems \cite{yu2018deep}, \cite{grohs2018proof}, \cite{sirignano2018dgm}. However, in this context, other well-known parametrizations are available \cite{bungartz2004sparse}. Besides that, current neural operators are inapplicable to high-dimensional problems.}

\begin{figure*}
    \centering
    \subfloat[]{
        \label{fig:bad_network}
        \includegraphics[width=0.45\textwidth]{pictures/bad_network-eps-converted-to.pdf}
    }
    \subfloat[]{
        \label{fig:inconsistent_FNO}
        \includegraphics[width=0.51\textwidth]{pictures/inconsistent_FNO_d-eps-converted-to.pdf}
    }
    \caption{(a) the output of neural network $N(x)$ computed on coarse and fine grids. On each subgrid, loss and gradients are zero, so the network provides the best (alas, pathological) approximation to $f(x) = 2x$ on the interval $[-1, 1]$ for all but the finest grid. (b) FNO trained to take the first derivative of random trigonometric polynomials on grid $2h$ is evaluated for the same functions on grid $h$. The difference between outputs for a particular function is on the right; statistics for $750$ functions is on the left. The average discrepancy is about $25\%$ of $L_2$ norm. Inconsistency of FNO results from aliasing of high harmonics produced by activation functions.}
    \label{fig:sno_motivation}
\end{figure*}

Second, the representation of the input function by sampling implicitly restrict the class of functions and possible operations with them. Indeed, it is known that for periodic functions sampling on grid with $2N+1$ points allows to resolve frequencies up to $N$ (Nyquist frequency).\footnote{For non-periodic functions the same argument can be made using Chebyshev polynomials.} All higher frequencies with $k>N$  are indistinguishable from lower frequencies with $k\leq N$ (this is known as aliasing \cite[Chapter 11, Definition 23]{boyd2001chebyshev}). For DeepONet this implies it is defined on trigonometric polynomials with degree $\leq N$. For FNO situation is more difficult, because even when the input function is fully resolved, nonlinear activations are going to introduce higher frequencies, so the lower frequencies will be distorted. As a result FNO $N(u(x), x)$ computed on grid $x_{2h}$ produces different result when computed on grid $x_{2h}$. This is illustrated in \cref{fig:inconsistent_FNO}, where the aliasing error $\left\|N(u(x_{2h}), x_{2h}) - \left[N(u(x_{h}), x_{h})\right]_{2h}\right\|$ ($\left[\cdot\right]_{2h}$ is a projection on grid $2h$) is computed.

Our solution to both of outlined problems is to detach the process of function approximation (interpolation) from the mapping between functional spaces, and explicitly fix the highest possible resolution. To do that, for both domain and codomain of neural operator we consider functions represented by finite series of the form $\sum c_n f_n(x)$ where $f_n$ are either Chebyshev or trigonometric polynomials. Now, all functions are equivalent to vectors of finite length, and we can realize mapping between functional spaces with (almost) vanilla feedforward neural network. Benefits of this strategy include: 1) reach set of operations (integration, differentiation, shift) available for Chebyshev and trigonometric polynomials; 2) lossless operations on functions; 3) interpolation on the fine mesh can be performed with arbitrary accuracy; 4) availability of global information about obtained functions (f.e. bounds on functions and their derivatives); 5) well-studied approximation properties; 6) good compression for smooth functions; 7) overall similarity with spectral and pseudospectral methods.

Main contributions are:
\begin{itemize}
    \item Analysis of aliasing error in FNO (\cref{subsection:Activation_functions_and_aliasing}, \cref{appendix:aliasing_and_relu}).
    \item Definition of super-resolution and generalization tests (\cref{subsection:Superresolution}, \cref{appendix:Additional_comments_on_super_resolution}).
    \item New type of Spectral Neural Operators (SNO), based on Chebyshev and Fourier series (\cref{section:Spectral_neural_operator}, \cref{appendix:details_on_training_and_architectures}).
    \item Comparison of SNO, DeepONet and FNO on $16$ problems (\cref{section:Numerical_examples}, \cref{appendix:benchmarks_for_additional_architectures}).
\end{itemize}

\section{Chebyshev and trigonometric polynomials}
\label{section:Chebyshev_and_trigonometric_polynomials}
As we stated, we are going to consider neural networks operating on functions of the form $f(x) = \sum c_n f_n(x)$, $x\in[-1, 1]$\footnote{As usual, basis functions  for higher dimensions are constructed as direct product of $D=1$ basis functions. From now on we consider $D=1$ having in mind that generalization to higher dimensions is straightforward.}, where $f_n(x)$ are either Chebyshev polynomials defined by reccurence \cite[Section 1.2.1]{mason2002chebyshev} $T_{n}(x) = 2x T_{n-1}(x) - T_{n-2}(x)$, $T_{0}(x)=1$, $T_{1}(x)=x$ or trigonometric polynomials $e^{i\pi kx}$, $k\in \mathbb{Z}$. This particular choice of basis warrants some explanation, which we provide in the present section.

First, coefficients $c_n$ can be efficiently approximated by sampling. The standard way to obtain coefficients is to use orthogonality of $f_n$ with respect to $L_2$ norm and compute the integral $\int f_{n}(x)f(x)dx$. For Chebyshev polynomials we can sidestep the integration using samples on Chebyshev grid $x_{k} = \cos(k\pi \big/ n),~k=0,\dots,n$, trigonometric identity $T_{n}(\cos(\theta)) = \cos(n\theta)$ and Discrete Cosine Transform  (DCT-I) \cite[Section 4.7]{mason2002chebyshev}. For trigonometric polynomials, coefficients can be obtained in a similar way but with sampling on the uniform grid and Fast Fourier Transform (FFT) \cite[Chapter 3]{trefethen2000spectral}. So the cost to recover $N$ coefficients is $O(N\log(N))$ plus the cost of $N$ evaluations of $f(x)$. The procedure also come with reliable adaptive algorithms which can chop the series based on accuracy requirements \cite{aurentz2017chopping}, \cite[Chapter 3]{boyd2014solving}, \cite{trefethen2007computing}. Moreover, the very same DCT-I and FFT can be used to perform interpolation both on Chebyshev and the uniform grid with the same complexity $O(N\log(N))$ \cite{boyd1992fast}. Altogether, we have a duality, between values on the grid and values of coefficients. This duality is significantly exploited in pseudospectral methods, when one needs to apply nonlinear transformation to given Chebyshev or trigonometric series.

The second reason is efficient compression of smooth functions. It is known that for functions with $m$ continuous derivatives $n$-th coefficient is $O(n^{-m})$ for both Chebyshev \cite[Theorem 5.14]{mason2002chebyshev} and Fourier series \cite[Chapter 2, Miscellanious theorems and examples]{zygmund2002trigonometric}.

Besides that, many operations such as integration, differentiation, shift (for Fourier series), multiplication of two series can be performed with machine accuracy \cite[Equation (2.49)]{mason2002chebyshev}, \cite[Equation (2.43)]{mason2002chebyshev}, \cite[Chapter 3]{trefethen2000spectral}, \cite{giorgi2011polynomial}.

All these properties combined make trigonometric and Chebyshev polynomials especially suitable for PDE solutions by spectral methods \cite{boyd2001chebyshev} and adaptive lossless computations with functions \cite{trefethen2007computing}. The later goal is fully realized in the Chebfun software.\footnote{\href{https://www.chebfun.org}{https://www.chebfun.org}} Chebfun demonstrates that computations with continuous functions can be performed with essentially the same guarantees as computations with floating-point numbers. So, in our opinion, Chebfun is a good role model for neural operators. In the present work, we use only a small part of Chebfun functionality, possible extensions are discussed in \cref{section:Conclusion_and_future_work}.

\section{Consequences of spectral representation}
\label{section:Consequences_of_spectral_representation}
Sufficiently smooth functions are indistinguishable from their Fourier and Chebyshev series if a reasonable number of terms is used. Given that, we restrict our attention to finite series and use duality between samples and coefficients to illustrate the accumulation of aliasing errors in FNO and discuss super-resolution. In what follows we mainly consider Fourier series and refer to $\phi(x; k) = \exp(i \pi k x)$ as harmonic with frequency $k$.

\subsection{Activation functions and aliasing}
\label{subsection:Activation_functions_and_aliasing}
Suppose our function $f(x)$ is (exactly) represented as Fourier series with $|k| < N$ terms. We can equivalently store values of the function on the uniform grid with $2N + 1$ points. However, when we apply activation function $\sigma(x)$, the resulting function $\sigma(f(x))$ may contain higher frequencies. If this is the case, composition $\sigma(f(x))$ can not be represented on the same grid with $2N+1$ points. The higher harmonics with $k>N$ can be computed on grid with $2N + 1$ points, but it is known that they are indistinguishable from lower harmonics with $k<N$ (see \cite[Chapter 11, Definition 23]{boyd2001chebyshev} for more details about Fourier aliasing and \cite[Chapter 4, Theorem 4.1]{trefethen2019approximation} for Chebyshev aliasing). If all operations are poinwise as in case of DeepONet or PiNN, it is not problematic. However, for FNO which uses global information, aliasing will present a problem. To quantify distortions introduced by activation function we define relative aliasing error $E_{a}$ for the Fourier series\footnote{For Chebyshev series norm is induced by scalar product with weight $1\big/\sqrt{1-x^2}$.}
\begin{equation}
    \label{aliasing_error}
    \sigma(f(x)) \equiv \sum_{i} c_i \phi_i,~E_a(f, \sigma)\equiv\frac{\left\|\sum_{|i|>N}c_i \phi_i\right\|_2}{\left\|\sigma(f(x))\right\|_2}.
\end{equation}
Aliasing error defined in \cref{aliasing_error} measures the norm of harmonics we cannot possibly resolve on the given grid relative to the norm of the function, transformed by activation. The following result gives aliasing error for rectifier and two extreme basis functions.

\begin{figure*}
    \centering
    \subfloat[]{
        \label{fig:FNO_aliasing}
        \includegraphics[width=0.48\textwidth]{pictures/oversampling_FNO_1-eps-converted-to.pdf}
    }
    \subfloat[]{
        \label{fig:super_resolution}
        \raisebox{-2pt}{\includegraphics[width=0.48\textwidth]{pictures/superresolution_error-eps-converted-to.pdf}}
    }
    \caption{(a) On the left, relative aliasing error for FNO $N(u, x)$, computed as $\left\|N(u(x_{2h}), x_{2h}) - \left[N(u(x_{h}), x_{h})\right]_{2h}\right\|$ normalized on $\left\|\left[N(u(x_{h}), x_{h})\right]_{2h}\right\|$, where $\left[\cdot\right]_{2h}$ is a projection on grid with spacing $2h$. FNO was trained on the grid with $n$ point ($x$ axis) on functions with frequencies $[0, 10]$. Aliasing error decreases when larger grids are used. On the right, relative test error for the input with $n$ point ($x$ axis) for FNO trained for input with $n=21$ points. Relative test errors is larger for finer grid, which means FNO has systematic bias caused by aliasing. (b) Neural operators were trained on input with harmonics $[k_{\min}, k_{\max}]$ and grid with $100$ points and evaluated on input with harmonics $[k_{\min} + \Delta k, k_{\max} + \Delta k]$ on grid with $300$ points. Relative test error sharply increases with $\Delta k$, super-resolution is not observed.}
    \label{fig:super_resolution_aliasing}
\end{figure*}

\begin{theorem}[Aliasing and ${\sf ReLU}$]
\label{theorem:aliasing_and_relu}
On the uniform grid with $2N + 1$ points and on Chebyshev grid with $N+1$ points for $E_a$ defined in \cref{aliasing_error} we have
\begin{equation}
    E_a(\cos(\pi N x), {\sf ReLU}) = E_a(T_N(x), {\sf ReLU}) = \frac{1}{\pi}\sqrt{\frac{\pi^2}{2} - 4}\simeq 0.31
\end{equation}
\end{theorem}
The proof can be found in \cref{appendix:aliasing_and_relu}. More results on aliasing for composition with smooth functions can be found in \cite{bergner2006spectral}. The aliasing error is quite substantial, but since all energy in the theorem above is confined in highest possible harmonics, in practice one can expect to have milder discrepancy.

\cref{theorem:aliasing_and_relu} predicts that in certain situations FNO has a systematic bias. Namely, if one trains FNO on an insufficiently fine grid, activation functions introduce distortions that FNO will learn to mitigate. When the grid is sufficiently refined, aliasing errors disappear, but since FNO was trained to mitigate them, it predicts output functions that differ from targets it was trained on. The precise discrepancy is hard to analyze but one can expect that for operators with smoother output aliasing will be milder. This hypothesis is confirmed by our experiments. Indeed, the right graph in \cref{fig:FNO_aliasing} shows that for the integration operator, which performs smoothing, aliasing leads to an approximately two-fold error increase, while for the differentiation operator the increase is five-fold. For both examples, FNO still produces good approximations, but one should be cautious using FNO for more complex problems, because in these applications FNO tends to be less accurate \cite{pathak2022fourcastnet}, \cite{grady2022towards}.

Our solution to the problem of aliasing errors is to decouple interpolation from the mapping between functional spaces. It is described in \cref{section:Spectral_neural_operator}. However, it is possible to decrease aliasing error even when FNO is used for interpolation. One simply needs to train (or fine-tune) on a sufficiently fine grid. The decrease of aliasing error in this scenario is illustrated in the left plot of \cref{fig:FNO_aliasing}.

\subsection{Super-resolution}
\label{subsection:Superresolution}

In a series of papers \cite{li2020fourier}, \cite{kovachki2021neural}, \cite{li2020neural} authors claim that neural operators can handle different resolution, and even demonstrate ``zero-shot super-resolution''. The definition of the later was not given, but ``zero-shot super-resolution'' is tested in two steps. First, FNO is trained on data on the uniform grid with $N$ points and reaches test error $E^{N}_{\text{test}}$. After that FNO is evaluated on the same data but on grid with $M>N$ points and demonstrates test error $E^{M}_{\text{test}} \simeq E^{N}_{\text{test}}$. This procedure is problematic for two reasons.

First, high-resolution data is used as an input during the second stage. It is abnormal, because in classical works on super-resolution several low-resolution images are combined to improve resolution \cite{tian2011survey} without the direct use of high-resolution data. Second, the presented procedure implies that resolution increases when the number of grid points increases. This is obviously not the case in general. To claim such a thing, one needs to specify the family of function she is working with. For example, consider functions given by trigonometric series with $2N+1$ terms. If we sample these functions on the uniform grid with $2N+1$ points, we can exactly recover all $2N+1$ coefficients and no further refinement of the grid is going to provide novel information about the function, i.e., to increase resolution.

Given that, we provide an improved super-resolution test that also measures the generalization properties of neural operators. At the first stage one trains neural operator on functions with harmonics $|k|\leq K$ on the grid with $N \geq 2K + 1$ points.\footnote{Here we refer to complex exponential $\phi_{k}(x)$ and uniform grid. The extension to the case of Chebyshev polynomials is obvious.} At the second stage, one evaluates neural operator on functions with some or all harmonics $|k| > K$ on the grid with $M \geq (2\max|k| + 1)$ points. Note, that our test does not take into account the properties of the operator to be learned. In some cases the test can produce meaningless results, see \cref{appendix:Additional_comments_on_super_resolution} for details.

\begin{wraptable}{r}{0.25\textwidth}
    \caption{Neural operators trained on high harmonics and evaluated on low harmonics.}
    \input{tables/low_frequencies}
    \label{table:low_frequencies}
\end{wraptable}
We perform this test on FNO, DeepONet and spectral neural operator (see \cref{section:Spectral_neural_operator}) trained to take derivative and perform transformation $f(x)\rightarrow f(x)f(x+1)$. The results and details of the experiments are in \cref{fig:super_resolution}.\footnote{DeepONet fails to learn the second operator, that is why data for DeepONet is absent in the second picture.} Clearly, no super-resolution is observed for neural operators. This is a frustrating experience because very simple integration and differentiation rules can pass this test effortlessly. The other consequence of experiment \cref{fig:super_resolution} is that all neural operators learn only frequencies that are present during training. This is true for low frequencies too. To illustrate this, we train FNO and spectral neural operator to take indefinite integral, using harmonics $k\in[15, 25]$ and evaluate both neural operators on data with lower harmonics $k\in[15-\Delta k, 25-\Delta k]$. The results are given in \cref{table:low_frequencies}. Again, we can see that neural operators are unable to perform well on low harmonics when they are not presented during the training stage.

\section{Spectral neural operator}
\label{section:Spectral_neural_operator}
Both results in \cref{table:low_frequencies} and in \cref{fig:super_resolution} suggest that neural operators perform mapping of finite-banded functions.\footnote{One may object that we considered only smooth data. This is, indeed, true, but as shown in \cite{de2022cost} both FNO and DeepONet produce Runge-like oscillations when applied to data with discontinuities, so with current versions of these operators non-smooth data is out of question.} Bands are fixed by the choice of architecture and training data. Then, the natural direction is to define neural operator with fixed number of harmonics (basis functions) for both input  and output function:
\begin{equation}
    \label{sno}
    \sum_{i}g_{i}(x) d_i = f_{\sf in}(x) \overset{\sf N}{\longrightarrow} f_{\sf out}(x) = \sum_{i}g_{i}(x) b_i,
\end{equation}
where $g_{i}(x)$ are either Chebyshev polynomial or complex exponential. So, neural network maps finite number of input coefficients $\left\{d_i\right\}$ to the finite number of output coefficients $\left\{b_i\right\}$. Using the duality between samples on the grid (Chebyshev or uniform) and values of coefficients we can alternatively map on-grid samples. We coin term {\em Spectral Neural Operator} (SNO) for any network that provide mapping specified in \cref{sno} and preserve the structure of the series.

It is important to understand \cref{sno} implies that SNO only realizes mapping between two functions given by truncated series. If one needs to compute these functions on the finer grid, Chebyshev or trigonometric interpolation should be use. In the same vein, Chebyshev/trigonometric smoothing (the chopping of the series) is used on the input function, if it is not conformable with a particular architecture of SNO.

\input{Tikz pictures/Chebyshev_LA}

The other important point is that series from \cref{sno} are truncated, so coefficients $d_i$ and $b_i$ can be stacked in finite-dimensional vectors. However, since these vectors represent sequences, linear algebra becomes slightly unusual. For example, it is possible to add two vectors even if their length does not coincide, or to multiply a vector by a matrix not conformable with it. The simplest way to understand this behavior is to think about sequences as vectors with an infinite number of trailing zeros that represent (zero) coefficients of truncated basis functions. They also can be considered as a simplified version of quasimatrices \cite{townsend2015continuous}. Standard linear algebra operations for sequences are illustrated in \cref{table:basic_operations_1}. It is also explained there, that, in the space of coefficients, linear operators quite appropriately becomes low-rank integral operators.

Except for unusual linear algebra, coefficients in \cref{sno} can be treated as ordinary vectors, so essentially any network that can map $\mathbb{C}^{n}$ to $\mathbb{C}^{m}$ qualifies as a neural operator.\footnote{The mapping can be between coefficients or, thanks to duality, between values of the function on the Chebyshev/uniform grid. When the later mapping is chosen, interpretation given in \cref{table:basic_operations_1} does not apply.} In the next section we gather benchmarks for SNO realized as ordinary feedforward neural networks, FNO and DeepONet.

\section{Numerical examples}
\label{section:Numerical_examples}

\begin{table*}[t]
    \centering
    \caption{Relative $L_2$ test errors for various problem. SNO(Ch) -- Spectral Neural Operator in Chebyshev basis, SNO(F) -- Spectral Neural Operator in the Fourier basis (in case of $2D$ problems, denoted by $xt$, Chebyshev basis is used along $t$ if the problem is not periodic), FNO \cite{li2020fourier} -- Fourier Neural Operator, DeepONet \cite{lu2019deeponet} -- Deep Operator Network.}
    \input{tables/truncated_big_table_with_all_results}
    \label{tab:big_table_with_all_results}
\end{table*}

Here we compare FNO \cite{li2020fourier}, DeepONet \cite{lu2019deeponet}, SNO(Ch) given by \cref{sno} in Chebyshev basis, SNO(F) given by \cref{sno} in Fourier basis. Precise architectures and details on training are described in \cref{appendix:details_on_training_and_architectures}. In addition, we perform benchmarks with more elaborate SNO architectures in \cref{appendix:benchmarks_for_additional_architectures}, and analyze sensitivity to initialization in \cref{appendix:sensitivity_to_initialization}. Code is available in \href{https://github.com/VLSF/SNO}{https://github.com/VLSF/SNO}.

In most of our benchmarks we used the following family of random functions
\begin{equation}
    \label{random_functions}
    \begin{split}
    &R(k_{\min}, k_{\max}, \sigma) = \left\{\frac{\mathcal{R}\left(\sum_{k=k_{\min}}^{k_{\max}}d_{k}^{j}e^{i\pi k x}\right)}{\sqrt{\sum_{k=k_{\min}}^{k_{\max}}\left|d_{k}^{j}\right|^{2}}}\right\}_{j},\,\mathcal{R}(d^{j}),\,\mathcal{I}(d^{j})\sim \mathcal{N}(0, \sigma^2),\,j\in[1, N]
    \end{split}
\end{equation}
where $\mathcal{R}(\cdot)$, $\mathcal{I}(\cdot)$ are real and imaginary part of complex number, $\mathcal{N}(a, b)$ is a normal distribution with mean $a$ and variance $b$. Operators that we used for training are described below.

\subsection{Operators}
\subsubsection{Integration, shift and differentiation}
\label{section:operators.subsection:integration_shift_and_differentiation}
First three test operators are
\begin{equation}
    \begin{split}
        &f(x) \rightarrow \int^{x} f(y) dy,\,f(y)\in R(1, 10, 2);\,f(x) \rightarrow f(x)f(x+1),\, f(x)\in R(0, 15, 2);\\
        &f(x) \rightarrow \frac{d}{dx}f(x),\, f(x)\in R(0, 10, k_{\max}),\,k_{\max}\in\left\{10, 20\right\}.
    \end{split}
\end{equation}

\subsubsection{Parametric ODE}
\label{section:operators.subsection:parametric_ode}
In this case we are interested in finding mapping $f(t)\rightarrow y(t)$, i.e., the solution to parametric family of ODEs (modification of DETEST problem $A4$ \cite{hull1972comparing})
\begin{equation}
    \frac{d y(t)}{dt} = y(t)f(t),\,y(0) = \exp\left(\int^{t=0} f(\tau) d\tau\right),\,f(t)\in R(1, 30, 2).
\end{equation}

\subsubsection{Elliptic equation}
\label{section:operators.subsection:elliptic_equation}
Here we considered two elliptic equations. In both cases, the goal is to predict the solution from the diffusion coefficient $k(x)$. In $1D$ we have
\begin{equation}
    \label{elliptic_1D}
    -\frac{d}{dx} k(x) \frac{d}{dx} u(x) = 1,\,u(-1)=u(1)=0,\,k(x) = 10(\tanh(g(x))+1) + 1,\,g(x)\in R(0, 20, 2).
\end{equation}
Diffusion coefficient changes from $1$ to $20$ but smoothly, so this equation is a variation of the benchmark (Darcy flow) used in \cite{li2020fourier}. In $2D$ we consider
\begin{equation}
    \label{elliptic_2D}
    \begin{split}
    &-\left(\frac{d}{dx} k(x)k(y) \frac{d}{dx} +\frac{d}{dy} k(x)k(y) \frac{d}{dy} \right)u(x, y) = 1,\,k(x) = \left(3(\tanh(g(x))+1) + 1\right),
    \end{split}
\end{equation}
with zero Dirichlet boundary conditions inside the square $[-1, 1]^2$ and $g(x),\,g(y)$ from $R(0, 20, 2)$.

\subsubsection{Burgers equation}
\label{section:operators.subsection:burgers_equation}
The goal is to infer solution of $ \frac{\partial}{\partial t} u(x, t) + \frac{1}{2}\frac{\partial}{\partial x}u(x, t)^2 - \nu \frac{\partial^2}{\partial x^2}u(x, t) = 0,\,u(x, 0) = f(x),\,x\in[-1,1]$, $f(x) \in R(0, 20, 2)$. In $1D$ we learn mapping $f(x)\rightarrow u(x, t=1)$, in $2D$ -- mapping $f(x)\rightarrow u(x, t),\,t\in[0, 1],\,x\in[-1, 1]$. We consider less viscous $\nu=0.01$ and more viscous $\nu=0.1$ scenarios.

\subsubsection{KdV equation}
\label{section:operators.subsection:kdv_equation}
In this benchmark we consider KdV equation
$\frac{\partial}{\partial t} u(x, t) + p\frac{\partial}{\partial x}u(x, t)^2 + \frac{\partial^3}{\partial x^3}u(x, t),\,u(x, 0) = f(x).$
More precisely, we use two exact solutions:
\begin{enumerate}
    \item Propagation of the single soliton ($p=1$) \cite{fornberg1978numerical}:
        \begin{equation}
            \label{soliton}
            u(x, t) = 3\left(\frac{a}{\cosh\left(a(x+x_0)\big/2 - a^3 t\big/2\right)}\right)^2,
        \end{equation}
    here $a$ defines amplitude and traveling speed and $x_0$ defines location of the maximum at $t=0$. We sample $a$ from uniform distribution on $[10, 25]$ and $x_0$ -- from uniform distribution on $[-1, 1]$. The goal is to predict $u(x, t=0.001)$ from $u(x, t=0)$ in $1D$ case, and $u(x, t),\,t\in[0,0.001]$ from $u(x, t=0)$ in $2D$ case.
    \item Interaction of two solitons ($p=6$) \cite{yoneyama1984korteweg}:
        \begin{equation}
            \label{two_solitons}
            \begin{split}
                &\phi_{i} = a_{i}(x+x_0^{i}) - 4a^3_{i} t,\\
                &\psi_1(x, t) = 2a_2^2(a_1^2 - a_2^2)\sinh\phi_1,\,
                \psi_2(x, t) = 2a_1^2(a_1^2 - a_2^2)\cosh\phi_2,\\
                &u(x, t) = \frac{\psi_1(x, t) + \psi_{2}(x, t)}{\left(a_1\cosh\phi_1\cosh\phi_2 - a_2\sinh\phi_1\sinh\phi_2\right)^2},\\
            \end{split}
        \end{equation}
        here $a_{i}$ defines the speed and amplitudes of solitons and $x^{i}_0$ their initial locations. We sample $a_1$ from uniform distribution on $[10, 25]$, $a_2=\chi a_1$, where $\chi$ is from uniform distribution on $[1/2, 1]$, the initial location was fixed to $x_0^{1}=-0.6$, $x_{0}^{2}=-0.5$. This way faster soliton is always behind the slower one and the neural operator needs to predict the whole interaction dynamics. For this problem time interval is $[0, 0.005]$, so the task is, again, to predict the solution in the final step, or on the whole spatio-temporal domain.
\end{enumerate}

\subsubsection{Nonlinear Schr\"odinger}
\label{section:operators.subsection:nonlinear_schrodinger}
Here we consider Kuznetsov-Ma breather \cite{kedziora2012second} given by
\begin{equation}
    \label{breather}
    \psi(x, t) = \left(\frac{-p^2\cos(\omega t) -2 i p\nu \sin(\omega t)}{2\cos(\omega t) - 2\nu \cosh(px\big/\sqrt{2})} - 1\right)e^{it},\,p=2\sqrt{\nu^2-1},\,\omega=p\nu.
\end{equation}
The goal is to predict $\left\|\psi(x, t)\right\|$ on $x\in[-1, 1]$, $t\in[0, 5]$ from initial conditions. The single parameter $\nu$ is sampled from uniform distribution on $[1.5, 3.5]$.

\subsection{Discussion of results}
\label{section:operators.subsection:discussion_of_results}
All results are gathered in \cref{tab:big_table_with_all_results}.

As we can see for many problems SNO performs better or comparable to FNO. The only exception is viscose Burgers equation with low viscosity. On this particular example both SNO networks failed to achieve reasonable test error. DeepONet perform much worse on our benchmarks. Especially, when high frequencies are present in the target. However note, that we use basic version of DeepONet. When more elaborate versions of DeepONet are used train error becomes comparable with the one of FNO \cite{lu2022comprehensive}.

\section{Conclusion and future work}
\label{section:Conclusion_and_future_work}
Using basic spectral techniques we describe the aliasing bias of FNO, discuss generalization properties and super-resolution of neural operators, and define SNO -- neural operator that operates on sequences of Chebyshev or Fourier coefficients. Benchmarks show that SNO is competitive when compared with other neural operators. So SNO appears to be a promising candidate for further investigations.

Limitations and natural lines for improvement are listed below.
\begin{enumerate}
    \item SNO in its current version is only applicable to smooth input and output data, because of the Gibbs phenomenon.\footnote{The same is true for DeepONet and FNO \cite{de2022cost}.} It is known that the Gibbs phenomenon can be completely resolved for spectral methods with a special reconstruction procedure \cite{gottlieb1992gibbs}. We plan to study this and related techniques in application to neural operators.
    \item Basis functions that we used are completely non-adaptive. In further research we plan to study how to implement an adaptive global basis using butterfly matrices \cite{dao2019learning}, which are known to be connected to fast special function transforms \cite{o2010algorithm}.
    \item The architecture that we used is very basic. It does not exploit in full the rich set of operations described in \cref{section:Chebyshev_and_trigonometric_polynomials}. We plan to implement and test substantially more complex SNO, possibly with adaptive operations defined in Chebfun\footnote{\href{https://www.chebfun.org}{https://www.chebfun.org}} \cite{trefethen2007computing}.
\end{enumerate}

\printbibliography[
heading=bibintoc,
title={Bibliography}
]


\section{Introduction}
In 1962, Kuhn gave us the idea that scientists adopt \emph{paradigms} (which are basically whole ways of doing science). Kuhn's model \cite{Kuhn1962} is that we have normal science governed by a single paradigm, where \emph{anomalies} build up, eventually leading to a \emph{crisis}, where scientists lose confidence in the paradigm. This leads to a \emph{revolution}, where one paradigm is replaced by another and normal science resumes. According to Kuhn, competing paradigms are \emph{incommensurable}, which basically means that they have \emph{no common measure for comparison}. This is one notion of incommensurability that we get from Kuhn. However, the concept of incommensurability goes beyond the absence of comparison bases. Some of the most interesting aspects of incommensurability were put forth by Paul Feyerabend (also in 1962). In fact, Feyerabend's notions of incommensurability \cite{feyerabend_1981} are not completely independent of Kuhn's notions, because he referred to Kuhn's then unpublished work while developing his ideas.

\section{Three Aspects of Incommensurability}

\noindent The concept of incommensurability (as developed by Kuhn and Feyerabend) can be broken down into three key ideas:
\begin{itemize}
    \item Methodological incommensurability.
    \item Observational incommensurability.
    \item Semantic incommensurability.
\end{itemize}

\noindent All these aspects are quite important and convey the essence of incommensurability. It therefore becomes necessary that we understand these three aspects before arguing for or against the role of incommensurability in physics\footnote{I will concern myself with the relevance of incommensurability only to physics.} (or any other science for that matter).

\subsection{Methodological Incommensurability}

The methodological aspect of incommensurability is exactly the one that I already mentioned. This is the claim that since each paradigm supplies its own standards of evaluation, there are no independent standards by which one can evaluate competing paradigms \cite{Hoyningen-Huene1990-HOYKCO}. In other words, one cannot step outside the realms of the competing paradigms to find an \emph{objective} standard of evaluation. Methodological incommensurability can be well understood by considering the case of gravitational theories before Einstein, i.e., Newton's and Descartes' theories of gravity.\\

\noindent One major problem with Newton's theory \cite{https://doi.org/10.48550/arxiv.1807.01791} was that it assumed \emph{non-locality}\footnote{Non-locality is the idea that information can instantaneously propagate from one point in space to another. This is a violation of the all-important postulates of relativity.} rather than explaining \emph{how} the gravitational force propagated. This \emph{action-at-a-distance} approach to gravity seemed quite mysterious to the physicists of Newton's time because they adopted the \emph{mechanistic approach to science}, where causes were supposed to be mechanistically linked to their effects.\\

\noindent A competing theory of gravity was offered by Descartes \cite{doi:10.1080/00033795900200038}, which claimed that all of space is filled with very fine matter, implying the non-existence of a \emph{truly empty space}. This matter tends to generate circular motion, which produces large \emph{vortices} in space. Descartes argued that the centrifugal force causes the fine matter to settle at the outer edges of these vortices, and subsequently, the inward pressure due to the fine matter pushes the heavier matter (like planets and stars) inwards. This pressure accounted for the effect of gravity.\\

\noindent Descartes' theory of vortices  is superior to that of Newton's in terms of providing a mechanistic explanation. On the other hand, Newton's theory produced \emph{better predictions}, which lead to its triumph over Descartes' theory. The interesting aspect of all this is that once Newton's theory was accepted into the scientific community, its lack of mechanistic explanations was no longer considered a problem. In fact, scientists altogether dropped the requirement that science must provide mechanistic explanations. Eventually, other scientists would have considered Newton's theory as an \emph{exemplar}\footnote{For instance, the theory of the electrostatic force developed by Coulomb is strikingly similar to Newton's theory of gravitation - both in terms of the equations, and in terms of the lack of a mechanistic explanation.} in the Kuhnian sense of the word.\\

\noindent The idea that we get from this historical account is that when the basis for evaluation of theories is changed, the problems of already existing theories tend to convert into virtues for any prospective theories. We see that due to changing methodological standards, Newton's theory became a paradigm - an exemplar against which one could judge other theories. The problems with Newton's theory were not solved by solving these problems but by altogether discarding the methodologies that produced those problems in the first place. Newton was more predictively accurate over Descartes while Descartes provided a mechanistic explanation. So, what is more important? Predictive accuracy or a mechanistic explanation? These are competing standards. One triumphed over the other, but did that have to be the case? It seems that scientists could have adopted for a different standard of evaluation and it is not clear that this choice would have been irrational or not.\\

\noindent In his later work, Kuhn allowed for certain standards to be widely accepted across paradigms. These include things like accuracy, consistency, scope, simplicity and fruitfulness for further research. But these virtues are imprecise and can be weighted in different ways. For instance, there is no algorithm to decide whether a theory is simple. Therefore, these virtues do not solve the problem of methodological incommensurability.

\subsection{Observational Incommensurability}

\noindent A rather more radical aspect of incommensurability is the observational aspect. This is more commonly known as the \emph{theory-ladenness of observation} \cite{Brewer2001-BRETTO}. The basic claim here is that the theories we accept alter the way we perceive the world. Standardly, we think of observations as providing a neutral arbiter between theories. For instance, if theory $A$ predicts $X$ and theory $B$ predicts $Y$, then it is straightforward to see which theory supports the observations. But according to Kuhn, the observations one makes are dependent on the theory that one accepts, which rejects the notion of the \emph{neutrality} of observations.\\

\noindent Norwood Russell Hanson gave a good example to understand this \cite{Hanson1958-HANPOD-2}. Suppose that Tycho Brahe and Johannes Kepler are watching the sun rise. Now, Brahe accepted the theory that the earth is the fixed center of the universe, while Kepler accepted that the earth orbits the sun. So, while observing the sun rise, Brahe would note down that the sun rises above a stationary horizon, while Kepler would note down that a moving horizon dips below the sun\footnote{One may argue that both these observations are correct, and that the only reason these observations differ is because the observers are in altogether different frames of reference (the relativistic argument), but note that Brahe and Kepler could be standing right next to each other (in the same frame of reference) and still make different observations.}. As Hanson puts it,

\begin{center}
   \emph{ There is more to seeing than meets the eyeball.} \cite{Hanson1958-HANPOD-2}
\end{center}

\noindent Kuhn provides another example for this in the context of particle physics:

\begin{center}
\noindent \emph{ Looking at a bubble-chamber photograph, the student\\ sees confused and broken lines, the physicist\\ a record of familiar subnuclear events.} \cite{Kuhn1962}
\end{center}

\noindent These examples, and more recent studies show that observations are indeed theory-laden. Hence, even if different paradigms had shared standards of evaluation, one would not be able to use evidence and logic alone to decide which one is better. Since scientists who adopt different paradigms see the world differently, the evidence is literally different for different paradigms. This implies that the decision to adopt to one of the competing paradigms is not a rational decision. In a nutshell, the data cannot settle the matter because the theories partially decide the data.\\

\noindent An obvious objection to this radical argument is that in order to say that there is no theory-neutral observation, one must show that this is true explicitly for all possible observations. These observations go well beyond what Kuhn and others have to offer. Other objections claim that the theory-ladenness of observations can be eliminated by careful reconsideration of the data. A seemingly good objection to observational incommensurability comes from Bogen and Woodward \cite{bogen1988saving}, who claim that theories are not tested against raw data, but against \emph{phenomena} which are unobservable and inferred from the data. Since phenomena are unobservable, they are not theory-laden, and hence one might get rid of observational incommensurability.

\subsection{Semantic Incommensurability}

\noindent The semantic aspect of incommensurability was popularized primarily by Feyerabend. The basic idea of semantic incommensurability is that the very meaning of the central terms in competing theories or paradigms are different (this is also called the problem of \emph{meaning invariance}), or as Feyerabend puts it, there is \emph{deductive disjointedness} between theories.\\

\noindent A great way to understand semantic incommensurability is by considering the development of physical ideas from Newton to Einstein. There is a clear disjointedness when it comes to the meaning of fundamental terms such as mass, energy, space, time, etc. in Newtonian mechanics as compared to Einstein's relativistic theories. In Newtonian mechanics, space and time are independent of each other and referred to as \emph{absolute space and time}. For Newton, space is flat and homogeneous, with a Euclidean structure imposed on it. Further, it is assumed that \emph{space and time exist of their own accord}, independent of mass and energy (matter). Einstein, on the other hand rejects the notion of absolute space and time, and puts forth the idea that space, time, mass and energy are all correlated \cite{Einstein:1905ve}. To quote John Wheeler,

\begin{center}
   \emph{Spacetime tells matter how to move; matter tells space-time how to curve}. \cite{thorne2000gravitation}
\end{center}

\noindent This is reflected in Einstein's famous field equations for gravitation \cite{Einstein:1915ca},

\begin{equation}
    R_{\mu\nu} - \frac{1}{2}Rg_{\mu\nu} + \Lambda g_{\mu\nu} = \frac{8\pi G}{c^4} T_{\mu\nu}
\end{equation}

\noindent where the term $T_{\mu\nu}$ on the RHS contains information about mass and energy, while everything on the LHS contains information about space and time - or rather spacetime. The equation mathematically proves Wheeler's statement. We see that Newtonian and Einsteinian ideas of fundamental entities have very little in common as the very meaning of terms change. Moreover, even the mathematical structure of these theories differ because Einstein imposes a non-Euclidean structure on spacetime, as opposed to Newton. We will discuss more examples of this across all of physics, when we consider physicists' arguments against incommensurability.

\section{What do Physicists have to say about Incommensurability?}

\noindent In 2002, the famous Russian trio of physicists - George Gamow, Dmitri Ivanenko and Lev Landau - published a paper titled \emph{World Constants and Limiting Transitions} \cite{article}. The paper talks about three things - \emph{theories} of physics (classical mechanics, quantum mechanics, etc.), \emph{world constants} (also called the fundamental constants of nature) and \emph{limiting transitions} (which is the idea that under different mathematical limits of these constants, one can make transitions from one theory to another).\\

\noindent Over the years, this idea has developed into a rich field of research which concerns itself with the unification of different physical theories into one all-encompassing framework called the theory of everything \cite{Manero2019-MANIOT-6}. The notion of world constants and limiting transitions has morphed into the pictorial representation called the cube of physics \cite{Padmanabhan2015}, which is depicted in figure 1. The theories in this model include classical mechanics (CM), quantum mechanics (QM), special relativity (SR), general relativity (GR), quantum field theory (QFT), Newtonian gravity (NG), non-relativistic quantum gravity (NRQG), and the theory of everything (TOE). These make up the vertices of the cube of physics. By themselves, these theories are just disjointed points having no correlation with each other. However, the effect of adding fundamental constants of nature to this model (i.e, the speed of light $c$, Planck's constant $\hbar$ and the universal gravitational constant $G$) is that of connecting these points. \\

\noindent In other words, we have theories, and we have bridges between these theories that are defined by the fundamental constants of nature. The way to think about this is as follows. Assume we have a consistent TOE which has all the constants ($c,\hbar,G$) present in its framework. Then in the absence of $\hbar$, the TOE reduces to GR, while in the absence of $c$, the TOE reduces to NRQG, and so on and so forth until we arrive at CM. If we run this process backwards, then we have a consistent framework to deduce a theory of everything from the already existing theories, i.e, one can unify all existing theories into a single theory.

\begin{figure}[!h]%
\centering
\includegraphics[scale=0.8]{coph.pdf}
\caption{The cube of physics describes physical theories and the transitions between these theories in terms of different conditions imposed on the fundamental constants of nature.}\label{fig1}
\end{figure}

\noindent The fact that the unification problem is probably pursued the most by physicists has something to say about the faith that physicists have on the cube of physics model. Physicists have shown that these transitions by \emph{the running of the constants} (as physicists like to call it) are true for already existing theories. For instance, SR does reduce to CM when $c \rightarrow \infty$. Hence, they inductively infer that the cube of physics is reliable, and that it will work for NRQG and TOE as well. What does this say about incommensurability? It seems as if it is possible to deductively infer one theory from another, unlike what Feyerabend said. At the outset, we get the idea that incommensurability may not be a problem after all. But after having a closer look, one can immediately find the flaws with this argument.\\

\noindent First of all, the conditions imposed on the fundamental constants are purely mathematical in nature, and carry no physical meaning. Consider the transition from quantum mechanics to classical mechanics. This is mathematically equivalent to the equation,

\begin{equation}\lim_{\hbar \to 0} \text{QM} = \text{CM}\end{equation}

\noindent This equation will only make sense to someone who does not care about the fundamental constants of nature. If one accepts that the Planck's constant $\hbar$ is a fundamental constant of nature, then the expression $\hbar \to 0$ makes no sense, because that would mean that the fundamental constant is no longer a constant, because its value changes! In this case, we conclude that such a limiting transition from quantum mechanics to classical mechanics is not possible. Hence, quantum mechanics and classical mechanics are incommensurable.\\

\noindent Secondly, even if we accept that $\hbar$ is not a fundamental constant of nature, and that its value could change, this does not prove anything about the truth of equation (3.1). It only implies that the equations of quantum mechanics (like Schrodinger's wave equation) reduce to the equations of classical mechanics (like Hamilton's equation). However, such mathematical deductions do not imply meaning invariance. The meanings of energy, mass, space, etc. in QM do not make a transition to their counterparts in CM in the limit where Planck's constant vanishes. So, even if $\hbar$ is not a constant, QM and CM are semantically incommensurable. In fact, this is exactly what physicists debate about in many conferences and lectures. In an ordinary undergraduate quantum mechanics lecture, no physicist will dare say that the concepts of QM have the same meaning as that of CM. On the other hand, they will (without hesitation) try to derive a mathematical relation between $[\hat{x},\hat{H}]$ (the quantum mechanical concept of motion) and \{$x,H$\} (the classical concept of motion).\\

\noindent These arguments will apply to any other transition described by the cube of physics, because they are all essentially of the same nature. But so far, we have only considered the physicists' case against semantic incommensurability. There are some popular arguments by physicists against methodological incommensurability as well. The most common one asserts the consideration of computational simplicity as a common measure to compare theories. Almost any classical physicist will say that the Hamiltonian and Lagrangian generalizations of mechanics are computationally much more simpler than the Newtonian counterpart, i.e, it is easier to solve problems using the action principle, rather than the force diagram approach of Newtonian mechanics. This is certainly true. In fact, the action principle has become a sort of consistency check for theoretical physics. Physicists say that if one cannot derive a theory from an action principle, then there is something wrong with the theory.\\

\noindent Again at the outset, this arguments discards methodological incommensurability. However, there is a flaw with this argument when one considers the subject of this argument. These arguments do not refer to different theories themselves, but to different mathematical formulations of the same theory. For instance, these arguments do not say that \emph{Newtonian mechanics} is better than \emph{general relativity} because of the computational simplicity of the former, but that the \emph{Lagrangian formulation} of classical mechanics (in the mathematical arena of configuration space) is better than the \emph{Newtonian formulation} (in the mathematical arena of physical space)\footnote{With quantum mechanics, the case is completely different, because its formulation in physical space (Bohmian mechanics) is altogether non-local.}. So, there is still no common base for comparison of theories. Methodological incommensurability persists.\\

\noindent Physicists often create a barrier of tests for new theories. This barrier includes concepts such as locality, Lorentz invariance, gauge invariance, etc. These are all certain mathematical tests for a theory. If the theory passes these tests, then it is a good theory. This might be used as an argument against methodological incommensurability, but the problem here is that this sounds exactly like the case with Newton's theory, where it became an exemplar for other theories. The new paradigms are tested against the laws of the existing paradigms, and thus there is no objective measure of comparison outside both these paradigms. Methodological incommensurability still persists.

\section{Concluding Remarks}

\noindent We have considered physicists' most powerful weapons against incommensurability and shown that these are either faulty weapons (i.e., they are inconsistent) or not weapons at all (i.e., they end up supporting the incommensurability thesis). These weapons are primarily used against methodological and semantic incommensurability. The validity of observational incommensurability is a bit unstable, but it can still be argued for in certain contexts.\\

\noindent The easiest way to argue for observational incommensurability is by comparing the notion of observations in classical and quantum physics. By definition, in classical physics, one is bound to observe classical objects (which have a certain minimum length scale or size). Likewise, in quantum physics, one is bound to observe subatomic particles or atoms (that are smaller than classical objects). This conveys the theory-ladenness of observations at the coarsest level. On the other hand, the very definition of a measurement or an observation is different in quantum mechanics \cite{RevModPhys.38.453,PhysRevD.2.2783}. An observation in quantum mechanics is supposed to permanently alter the state of the system that one is observing. This weirdness is not to be found in classical physics.\\

\noindent This might seem like a good argument for observational incommensurability, and one can always dig a little deeper by looking at more specific examples. But, as aforementioned, the only way to argue for observational incommensurability (at least to the best of my knowledge) is to prove it explicitly for every possible scenario. The answer to the question posed in the title of this paper is that even though methodological and semantic incommensurability can be considered to be existing \emph{globally} across all of physics, one is only allowed to argue for the \emph{local} existence of observational incommensurability within a few physical models, unless its existence is proved across all possible models.

\bibliographystyle{IEEEtran}

\section{Introduction }
Oxford University Press has developed this authoring template to help authors in preparing articles. Authors are encouraged to use this template to produce \LaTeX\ manuscripts which conform to IMA Journal styles. This document is a manual for authors to help using this template during article preparation. It features general guidelines and contains descriptions regarding various elements that can be used while preparing manuscripts. Authors are requested to refer to the file ``README'' for details of files available for reference.

Please utilize the IMA authoring template to the maximum, rather than
adding further packages or macros if possible. Please use semantic mark-up instead of formatting markup where there is a choice, for instance using \verb+\emph+ rather than \verb+\it+ or \verb+\textit+ to indicate emphasis. This will ease the production process and avoids misinterpretation of your LATEX code.

This documentation is not intended to give an introduction to \LaTeX. For questions concerning \TeX\ systems/installations or the \LaTeX\ mark-up language in general, please visit \url{http://tug.ctan.org/} or any other \TeX\ user group worldwide. The essential reference for \LaTeX\ is Mittelbach F., Goossens M. (2004) \textit{The \LaTeX\ Companion. 2nd edn.}, but there are many other good books about \LaTeX.

\section{Quick reference guide}
\begin{itemize}
\item How to download and install - \url{https://miktex.org/download}
\item Any system requirements - \url{https://miktex.org/kb/prerequisites-2-9}
\item Quick steps to get started -\newline \url{http://users.dickinson.edu/~richesod/latex/latexcheatsheet.pdf}\newline
\url{http://wch.github.io/latexsheet/latexsheet.pdf}\newline
\url{https://www.overleaf.com/learn}
\end{itemize}

\section{How to start and prepare your article  }
It is assumed that you possess basic knowledge of  \LaTeX. Unless using \LaTeX\ online (for example via Overleaf.com) ensure that you have \LaTeX2e version installed on your computer. You are provided with a class file ``ima-authoring-template.cls''. This template can be kept with your manuscript files. Note that the class file depends on the following packages which are standard and are available along with \LaTeX installation:

\begin{center}\begin{tabular}{|l|l|l|l|l|}\hline
graphicx & multirow & amsmath & amssymb & amsfonts \\\hline
array & flushend & stfloats & color & xcolor \\\hline
rotating & chngpage & totcount & fix-cm & algorithm \\\hline
algorithmicx & algpseudocode & listings & url & crop \\\hline
\end{tabular}\end{center}

Apart from the above-listed packages, the following additional packages are used in the class file for providing add-on functionalities to the template:

\begin{center}\begin{tabular}{|c|c|c|c|c|}\hline
subfloat & subfig & appendix & hyperref & footnote \\\hline
mathrsfs & natbib & wrapfig & amsthm & apacite \\\hline
\end{tabular}\end{center}

To learn more about the underlying packages, please read the respective documentations (try, e.g., texdoc [package name] at your shell prompt or visit \url{http://tug.ctan.org/}).

\medskip

\noindent You are given with sample tex file: ima-authoring-template.tex

\medskip

Please save a copy of the template file before you start editing the file as per your requirement. These samples contain the lines for calling class files, the preamble area, and the start/end of document where major sample elements required for an article are placed. Comments are included for each element and are self-explanatory. You can add your actual manuscript content in place of these sample elements. The standard structure of each element for an article is explained in detail in the following sections.

To use IMA authoring template, put all the package files in your working directory, edit the file ``sample template file'' in your preferred text editor, and run \LaTeX\ as usual. The resulting layout is similar but not identical to the layout of the final article.

Please note that you are not responsible for any final page layout. It is not necessary (and is even sometimes undesirable) to do any fine-tuning with commands like \verb+\break+, \verb+\pagebreak+, \verb+\vspace+, \verb+\clearpage+, etc. Please use semantic mark-up as far as possible and avoid additional formatting commands.

\section{Package features and some important settings}

\subsection{Language}
English is the default language used for typesetting rules.

\subsection{Fonts}
Please refrain from using custom fonts.

\medskip

\noindent Text fonts: Unlike the final published version, the authoring template uses non-commercial fonts: Times. These fonts are free version of the PostScript standard fonts and are supplied as part of all standard \TeX\ distributions.

\medskip

\noindent Math fonts: The standard Computer Modern math fonts are used.

\section{Preamble}
The preamble part comes between the document class line---\newline\verb+\documentclass{...}+---and the beginning of your document---\verb+\begin{document}+. Use this preamble area to include additional packages and customized macros if any.

\section{Details on document class options}

\begin{enumerate}
\item[a.] Options available to select below IMA journals:
\begin{itemize}
\item IMAIAI
\item IMAMAN
\item IMAMAT
\item IMAMCI
\item IMAMMB
\item IMANUM
\item IMATRM
\end{itemize}
\item[c.] namedate --- for authoryear citation style; default reference citation style is numbered reference style
\item[d.] webpdf --- for cropped paper size in the PDF output
\item[e.] unnumsec --- to get unnumbered section heads
\end{enumerate}

\section{Major structures/elements}

Article contents are divided into three main elements--front matter, main matter, and back matter. The elements preceding the \verb+\maketitle+ tag are considered as front matter elements and the elements placed below the \verb+\maketitle+ tag are main matter elements. Elements found after the section ``Conclusion'' are considered as back matter elements.

\begin{center}\tabcolsep2pt\fontsize{6}{9}\selectfont\begin{tabular}{|l|l|l|}\hline
\textbf{Front matter}&   \textbf{Main matter}&    \textbf{Back matter}\\
\hline
\verb+\title{...}+&     \verb+\section{...}+&   \verb+\begin{appendices}...\end{appendices}+\\\hline
\verb+\author{...}+&    \verb+\subsection{...}+&        \verb+\bibliographystyle{...}+\\\hline
\verb+\authormark{...}+&  \verb+\subsubsection{...}+&     \verb+


\section{0pt}{12pt plus 4pt minus 2pt}{0pt plus 2pt minus 2pt}
\usepackage[a4paper, total={7.0in, 10.5in}]{geometry}

\title{Pushing compute and AI onto detector silicon\\
AI@DOE Roundtable Whitepaper 2021}

\author[1]{Antonino Miceli \thanks{amiceli@anl.gov}}
\author[2]{Kazutomo Yoshii \thanks{kazutomo@mcs.anl.gov}}
\author[3]{Ian T. Foster \thanks{foster@anl.gov}}
\affil[1]{X-ray Science Division,\
Argonne National Laboratory, Lemont, IL USA}
\affil[2]{Mathematics \& Computer Science Division,\
Argonne National Laboratory, Lemont, IL USA}
\affil[3]{Data Science \& Learning Division,\
Argonne National Laboratory, Lemont, IL USA}
\date{12/08/2021}

\begin{document}
\maketitle
\vspace{-30pt}

\section{Challenge}
Pixel detectors are at the heart of advances in imaging. Astronomy has been revolutionized by the use of high-throughput detectors for surveying supernovas, shedding light on expansion rates of the universe, and exoplanet searches using the wobble of faint stars. Biochemistry has been transformed first by developing megapixel detectors for X-ray crystallography, and more recently, by using high frame rate detectors to correct for drift and enable atomic resolution structure determination in cryo-electron microscopy. Materials scientists are gaining new insights into subtle electron bonding arrangements using high dynamic range electron detectors as part of phasing coherent diffraction patterns.  Beyond scientific advances, the photography industry has been revolutionized by the transition to electronic image detectors.

However, in order to take full advantage of the DOE’s billion-dollar investments into the next-generation research infrastructure (e.g., exascale, light sources, colliders), advances are required not only in detector technology but also in computing and specifically AI. Let us consider an example from X-ray science. Nanoscale X-ray imaging is a crucial tool to enable a wide range of scientific explorations from materials science and biology to mechanical and civil engineering. The next-generation light sources will increase the X-ray beam brightness and coherent flux by 100 to 1,000 times. In order to image larger samples, the continuous frame rate of pixel array detectors must be increased, approaching 1 MHz, which requires several Tbps (aggregated) to transfer pixel data out to a data acquisition system. Using 65-nm CMOS technology, an optimistic raw data rate off such a chip is about 100-200 Gbps. However, a continuous 1 MHz detector with only $256 \times 256$  pixels at 16-bit resolution, for example, will require 1,000 Gbps (i.e., 1 Tbps) bandwidth off the chip!  It is impractical to have multiple high-speed transceivers running in parallel to provide such bandwidth and represents the first data bottleneck. New approaches are necessary to reduce the data size by performing data compression or AI-based feature extraction directly inside a detector silicon chip in a streaming manner before sending it off-chip.

The dividing line between edge and exascale computing is often ambiguous. However for scientific imaging detectors, the ultimate location to incorporate edge compute is directly on the front-end detector silicon.  Many of these scientific detectors are constructed from application-specific integrated circuits (ASICs). Until recently, pixel detector ASICs have been used mainly for analog signal processing of the charge from the sensor layer and the transmission of raw pixel data off the detector ASIC. However, with the availability of more advanced ASIC technology nodes for scientific application, more digital functionality from the computing domains (e.g., compression and feature extraction) can be  integrated directly into the detector ASIC (e.g., on the chip's periphery or edge) to overcome off-chip bandwidth limitation and ultimately increase data velocity. This increase in data velocity will not only enable observation phenomena at faster timescales, but also enable smarter experiments to more quickly focus on important regions of interest and detect rare events. We should stress that an increase in data velocity can only be achieved by embedding AI into the detector silicon itself. External AI accelerators (e.g., GPUs, TPUs, etc) can only help with downstream bottlenecks and will not fundamentally increase data velocity.

\section{Opportunity}
The scientific pixel detectors are typically built using ASICs. In general, integrated circuits are at the core of smartphones and cameras, the internet of things, computers, cloud data centers, and a new wave of artificial intelligence processors. What has enabled the proliferation of integrated circuits beyond traditional computing units (i.e., CPUs) is the “pure-play” semiconductor foundry business model, where even small-scale companies can develop designs that are then fabricated using billion-dollar semiconductor fabrication facilities. This resulted in a new opportunity for fabless semiconductor companies to emerge which focus only on integrated circuit design. The academic research community is fortunate to have access to these foundries through third-party brokers who facilitate cost-effective prototyping of small designs and full-scale fabrications at  advanced technology nodes. With the availability of these advanced technology nodes, the number of digital logic resources available increases. The use of digital circuitry also opens up the possibility of incorporating  digital data processing (e.g., compression, machine learning) directly on the detector ASIC.

The current generation of electron and X-ray pixel detectors and those on the horizon work at kHz frame rates in continuous mode. Pixel detectors also exist that capture images at $\sim$~10 MHz frame rates in burst mode. For example, the EU-XFEL’s AGIPD and LPD detectors and Cornell University’s Keck-PAD use on-chip analog storage to capture bursts of 8 to 384 images at MHz rates but  are subsequently read out slowly.  These detector projects demonstrate that the current generation of analog front ends can handle $\sim$~10 MHz frame rates. The common bottleneck of these detectors is the limited data bandwidth off the detector ASIC resulting from the analog design flow. Digital technology can increase the off-chip bandwidth if we follow a more modern mixed-signal ASIC design flow methodology using sub-100 nm CMOS processes.

While recent exploratory work to develop compute on detector ASICs is promising~\cite{Hammer_2021,Yoshii_codesign,Strempfer:2021uo}, much more is required. A co-design process that brings together ASIC designers, computer and application scientists is critical to make significant progress. In addition, a DOE-wide effort is needed to merge compute with detector ASICs. However, a big hurdle to cross-lab efforts are the significant legal issues with NDAs from chip foundries as well as the very high entry to barrier of commercial electronic design automation (EDA) tools since academic pricing is not offered to DOE labs. Ideally, one would like to have affordable EDA tools and NDAs where all the DOE labs can share circuit designs which each other, just like computer scientists have done with open-source software. A substantial coordinated effort is required to address these issues. Technological innovation should be centered around ultrafast inference architecture designs and software technologies to assist hardware development for detectors. This research effort on ultrafast inference architecture could also benefit future high-performance computation as workloads for processing units (e.g., FPU) shift from conventional numerical computing to AI workloads.

\section{Timeliness or maturity}
We believe that now is the perfect time to invest in the development of a new generation of AI-accelerated detectors. Billions of dollars are being spent on next-generation light source and collider facilities in the US (ALS-U, APS-U, CHESS-U, LCLS-II, EIC) and around the world (MAX-IV, ESRF-U, PETRA-IV, Beijing Light Source, etc.) as well as at leadership computing facilities. While the brightness increases with these new  accelerator upgrades will be significant, these new sources will not reach their full potential if limited to using detectors of the type available today, or those presently seen on the horizon. In addition, application scientists have increasingly turned to AI to analyze data~\cite{sanaullah2018real,Abeykoon2019,Zhengchun2019,Cherukara:2020il,liu2020braggnn}. AI-accelerated workflows have been shown not only to be fast enough to keep up with experiments, but also to overcome experimental restrictions of conventional methods. Now is the time to merge compute, AI and detectors!

\bibliographystyle{ieeetr}