Operator Learning Meets Numerical Analysis: Improving Neural Networks
  through Iterative Methods
Deep neural networks, despite their success in numerous applications, often
function without established theoretical foundations. In this paper, we bridge
this gap by drawing parallels between deep learning and classical numerical
analysis. By framing neural networks as operators with fixed points
representing desired solutions, we develop a theoretical framework grounded in
iterative methods for operator equations. Under defined conditions, we present
convergence proofs based on fixed point theory. We demonstrate that popular
architectures, such as diffusion models and AlphaFold, inherently employ
iterative operator learning. Empirical assessments highlight that performing
iterations through network operators improves performance. We also introduce an
iterative graph neural network, PIGN, that further demonstrates benefits of
iterations. Our work aims to enhance the understanding of deep learning by
merging insights from numerical analysis, potentially guiding the design of
future networks with clearer theoretical underpinnings and improved
performance.
Emanuele Zappala, Daniel Levine, Sizhuang He, Syed Rizvi, Sacha Levy, David van Dijk
# Operator Learning Meets Numerical Analysis: Improving Neural Networks through Iterative Methods

###### Abstract

Deep neural networks, despite their success in numerous applications, often function without established theoretical foundations. In this paper, we bridge this gap by drawing parallels between deep learning and classical numerical analysis. By framing neural networks as operators with fixed points representing desired solutions, we develop a theoretical framework grounded in iterative methods for operator equations. Under defined conditions, we present convergence proofs based on fixed point theory. We demonstrate that popular architectures, such as diffusion models and AlphaFold, inherently employ iterative operator learning. Empirical assessments highlight that performing iterations through network operators improves performance. We also introduce an iterative graph neural network, PIGN, that further demonstrates benefits of iterations. Our work aims to enhance the understanding of deep learning by merging insights from numerical analysis, potentially guiding the design of future networks with clearer theoretical underpinnings and improved performance.

## 1 Introduction

Deep neural networks have become essential tools in domains such as computer vision, natural language processing, and physical system simulations, consistently delivering impressive empirical results. However, a deeper theoretical understanding of these networks remains an open challenge. This study seeks to bridge this gap by examining the connections between deep learning and classical numerical analysis.

By interpreting neural networks as operators that transform input functions to output functions, discretized on some grid, we establish parallels with numerical methods designed for operator equations. This approach facilitates a new iterative learning framework for neural networks, inspired by established techniques like the Picard iteration.

Our findings indicate that certain prominent architectures, including diffusion models, AlphaFold, and Graph Neural Networks (GNNs), inherently utilize iterative operator learning (see Figure 1). Empirical evaluations show that adopting a more explicit iterative approach in these models can enhance performance. Building on this, we introduce the Picard Iterative Graph Neural Network (PIGN), an iterative GNN model, demonstrating its effectiveness in node classification tasks.

In summary, our work:
* Explores the relationship between deep learning and numerical analysis from an operator perspective.
* Introduces an iterative learning framework for neural networks, supported by theoretical convergence proofs.
* Evaluates the advantages of explicit iterations in widely-used models.
* Presents PIGN and its performance metrics in relevant tasks.
* Provides insights that may inform the design of future neural networks with a stronger theoretical foundation.

The remainder of this manuscript is organized as follows: We begin by delving into the background and related work to provide the foundational understanding for our contributions. This is followed by an introduction to our theoretical framework for neural operator learning. Subsequently, we delve into a theoretical exploration of how various prominent deep learning frameworks undertake operator learning. We conclude with empirical results underscoring the advantages of our proposed framework.

## 2 Background and Related Work

Numerical Analysis.Numerical analysis is rich with algorithms designed for approximating solutions to mathematical problems. Among these, the Banach-Caccioppoli theorem is notable, used for iteratively solving operator equations in Banach spaces. The iterations, often called Fixed Point iterations, or Picard iterations, allow to solve an operator equation approximately, in an iterative manner. Given an operator \(T\), this approach seeks a function \(u\) such that \(T(u)=u\), called a fixed point, starting with an initial guess and refining it iteratively.

The use of iterative methods has a long history in numerical analysis for approximate solutions of intractable equations, for instance involving nonlinear operators. For example, integral

Figure 1: Overview of iterative framework. (A) Popular architectures which incorporate iterative components in their framework. (B) Convergence behavior of an iterative solver. (C) Behavior of iterative solver converging to a fixed point in the data manifold.
equations, e.g. of Urysohn and Hammerstein type, arise frequently in physics and engineering applications and their study has long been treated as a fixed point problem [13, 14, 15].

Convergence to fixed points can be guaranteed under contractivity assumptions by the Banach-Caccioppoli fixed point theorem [1]. Iterative solvers have also been crucial for partial differential equations and many other operator equations [16].

Operator Learning.Operator learning is a class of deep learning methods where the objective of optimization is to learn an operator between function spaces. Examples and an extended literature can be found in [17, 18]. The interest of such an approach, is that mapping functions to functions we can model dynamics datasets, and leverage the theory of operators. When the operator learned is defined through an equation, e.g. an integral equation as in [19], along with the training procedure we also need a way of solving said equation, i.e. we need a solver. For highly nonlinear problems, when deep learning is not involved, these solvers often utilize some iterative procedure as in [16]. Our approach here brings the generality of iterative approaches into deep learning by allowing to learn operators between function spaces through iterative procedure used in solving nonlinear operator equations.

Transformers.Transformers ([23, 17, 18]), originally proposed for natural language processing tasks, have recently achieved state-of-the-art results in a variety of computer vision applications ([19, 18, 16, 17, 20, 2]). Their self-attention mechanisms make them well-suited for tasks beyond just sequence modeling. Notably, transformers have been applied in an iterative manner in some contexts, such as the "recycling" technique used in AlphaFold2 [15].

AlphaFold.DeepMind's AlphaFold [19] is a protein structure prediction model, which was significantly improved in [15] with the introduction of AlphaFold2 and further extended to protein complex modeling in AlphaFold-Multimer [14]. AlphaFold2 employs an iterative refinement technique called "recycling", which recycles the predicted structure through its entire network. The number of iterations was increased from 3 to 20 in AF2Complex [1], where improvement was observed. An analysis of DockQ scores with increased iterations can be found in [1]. We only look at monomer targets, where DockQ scores do not apply and focus on global distance test (GDT) scores and root-mean-square deviation (RMSD).

Diffusion Models.Diffusion models were first introduced in [13] and were shown to have strong generative capabilities in [19] and [20]. They are motivated by diffusion processes in non-equilibrium thermodynamics [14] related to Langevin dynamics and the corresponding Kolmogorov forward and backward equations. Their connection to stochastic differential equations and numerical solvers is highlighted in [13, 14, 15, 16]. We focus on the performance of diffusion models at different amounts of timesteps used during training, including an analysis of FID [12] scores.

Graph Neural Networks (GNNs).GNNs are designed to process graph-structured data through iterative mechanisms. Through a process called message passing, they repeatedly aggregate and update node information, refining their representations. The iterative nature of GNNs was explored in [12], where the method combined repeated applications of the same GNN layer using confidence scores. Although this shares similarities with iterative techniques, our method distinctly leverages fixed-point theory, offering specific guarantees and enhanced performance, as detailed in Section 5.1.

## 3 Iterative Methods for Solving Operator Equations

In the realm of deep learning and neural network models, direct solutions to operator equations often become computationally intractable. This section offers a perspective that is applicable to
machine learning, emphasizing the promise of iterative methods for addressing such challenges in operator learning. We particularly focus on how the iterative numerical methods converge and their application to neural network operator learning. These results will be used in the Appendix to derive theoretical convergence guarantees for iterations on GNNs and Transformer architectures, see Appendix A.

### Setting and Problem Statement

Consider a Banach space \(X\). Let \(T:X\longrightarrow X\) be a continuous operator. Our goal is to find solutions to the following equation:

\[\lambda T(x)+f=x, \tag{1}\]

where \(f\in X\) and \(\lambda\in\mathbb{R}-\{0\}\) is a nontrivial scalar. A solution to this equation is a fixed point \(x^{*}\) for the operator \(P=\lambda T+f\):

\[\lambda T(x^{*})+f=x^{*}. \tag{2}\]

### Iterative Techniques

It is clear that for arbitrary nonlinear operators, solving Equation (1) is not feasible. Iterative techniques such as Picard or Newton-Kantorovich iterations become pivotal. These iterations utilize a function \(g\) and progress as:

\[x_{n+1}=g(T,x_{n}). \tag{3}\]

Central to our discussion is the interplay between iterative techniques and neural network operator learning. We highlight the major contribution of this work: By using network operators iteratively during training, convergence to network fixed points can be ensured. This approach uniquely relates deep learning with classical numerical techniques.

### Convergence of Iterations and Their Application

A particular case of great interest is when the operator \(T\) takes an integral form and \(X\) represents a function space, our framework captures the essence of an integral equation (IE). By introducing \(P_{\lambda}(x)=\lambda T(x)+f\), we can rephrase our problem as a search for fixed points.

We now consider the problem of approximating a fixed point of a nonlinear operator. The results of this section are applied to various deep learning settings in Appendix A to obtain theoretical guarantees for the iterative approaches.

**Theorem 1**.: _Let \(\epsilon>0\) be fixed, and suppose that \(T\) is Lipschitz with constant \(k\). Then, for all \(\lambda\) such that \(|\lambda k|<1\), we can find \(y\in X\) such that \(||\lambda T(y)+f-y||<\epsilon\) for any choice of \(\lambda\), independently of the choice of \(f\)._

Proof.: Let us set \(y_{0}:=f\) and \(y_{n+1}=f+\lambda T(y_{n})\) and consider the term \(||y_{1}-y_{0}||\). We have

\[||y_{1}-y_{0}||=||\lambda T(y_{0})||=|\lambda|||T(y_{0})||.\]

For an arbitrary \(n>1\) we have

\[||y_{n+1}-y_{n}||=||\lambda T(y_{n})-\lambda T(y_{n-1})||\leq k|\lambda|||y_{ n}-y_{n-1}||.\]

Therefore, applying the same procedure to \(y_{n}-y_{n-1}=T(y_{n-1})-T(y_{n-2})\) until we reach \(y_{1}-y_{0}\), we obtain the inequality

\[||y_{n+1}-y_{n}||\leq|\lambda|^{n}k^{n}||T(y_{0})||.\]

Since \(|\lambda|k<1\), the term \(|\lambda|^{n}k^{n}||T(y_{0})||\) is eventually smaller than \(\epsilon\), for all \(n\geq\nu\) for some choice of \(\nu\). Defining \(y:=y_{\nu}\) for such \(\nu\) gives the following

\[||\lambda T(y_{\nu})+f-y_{\nu}||=||y_{\nu+1}-y_{\nu}||<\epsilon.\]
The following now follows easily.

**Corollary 1**.: _Consider the same hypotheses as above. Then Equation 1 admits a solution for any choice of \(\lambda\) such that \(|\lambda|k<1\)._

Proof.: From the proof of Theorem 1 it follows that the sequence \(y_{n}\) is a Cauchy sequence. Since \(X\) is Banach, then \(y_{n}\) converges to \(y\in X\). By continuity of \(T\), \(y\) is a solution to Equation 1. 

Recall that for nonlinear operators, continuity and boundedness are not equivalent conditions.

**Corollary 2**.: _If in the same situation above \(T\) is also bounded, then the choice of \(\nu\) of the iteration can be chosen uniformly with respect to \(f\), for a fixed choice of \(\lambda\)._

Proof.: From the proof of Theorem 1, we have that

\[||y_{n+1}-y_{n}||\leq|\lambda|^{n}k^{n}||T(y_{0})||=|\lambda|^{n}k^{n}||T(f)||.\]

If \(T\) is bounded by \(M\), then the previous inequality is independent of the element \(f\in X\). Let us choose \(\nu\) such that \(|\lambda|^{n}k^{n}<\epsilon/M\). Then, suppose \(f\) is an arbitrary element of \(X\). Initializing \(y_{0}=f\), \(y_{\nu}\) will satisfy \(||\lambda T(y_{\nu})+f-y_{\nu}||<\epsilon\), for any given choice of \(\epsilon\). 

The following result is classic, and its proof can be found in several sources. See for instance Chapter 5 in [1].

**Theorem 2**.: _(Banach-Caccioppoli fixed point theorem) Let \(X\) be a Banach space, and let \(T:X\longrightarrow X\) be contractive mapping with contractivity constant \(0<k<1\). Then, \(T\) has a unique fixed point, i.e. the equation \(T(x)=x\) has a unique solution \(u\) in \(X\). Moreover, for any choice of \(u_{0}\), \(u_{n}=T^{n}(u_{0})\) converges to the solution with rate of convergence_

\[||u_{n}-u|| < \frac{k^{n}}{1-k}||u_{0}-u_{1}||, \tag{4}\] \[||u_{n}-u|| < \frac{k}{1-k}||u_{n-1}-u_{n}||. \tag{5}\]

The possibility of solving Equation 1 with different choices of \(f\) is particularly important in the applications that we intend to consider, as it is interpreted as the initialization of the model. While various models employ iterative procedures for operator learning tasks implicitly, they lack a general theoretical perspective that justifies their approach. Several other models can be modified using iterative approaches to produce better performance with lower number of parameters. We will give experimental results in this regard to validate the practical benefit of our theoretical framework.

While the iterations considered so far have a fixed procedure which is identical per iteration, more general iterative procedures where the step changes between iterations are also diffused, and this can be done also adaptively.

### Applications

Significance and Implications.Our results underscore the existence of a solution for Equation 1 under certain conditions. Moreover, when the operator \(T\) is bounded, our iterative method showcases uniform convergence. It follows that ensuring that the operators approximated by deep neural network architectures are contractive, we can introduce an iterative procedure that will allow us to converge to the fixed point as in Equation 2.

Iterative Methods in Modern Deep Learning.In contemporary deep learning architectures, especially those like Transformers, Stable Diffusion, AlphaFold, and Neural Integral Equations, the importance of operator learning is growing. However, these models, despite employing iterative techniques, often lack the foundational theoretical perspective that our framework provides. We will subsequently present experimental results that vouch for the efficacy and practical advantages of our theoretical insights.
Beyond Basic Iterations.While we have discussed iterations with fixed procedures, it is imperative to highlight that more general iterative procedures exist, and they can adapt dynamically. Further, there exist methods to enhance the rate of convergence of iterative procedures, and our framework is compatible with them.

## 4 Neural Network Architectures as Iterative Operator Equations

In this section, we explore how various popular neural network architectures align with the framework of iterative operator learning. By emphasizing this operator-centric view, we unveil new avenues for model enhancements. Notably, shifting from implicit to explicit iterations can enhance model efficacy, i.e. through shared parameters across layers. A detailed discussion of the various methodologies given in this section is reported in Appendix B. In the appendix, we investigate architectures such as neural integral equations, transformers, AlphaFold for protein structure prediction, diffusion models, graph neural networks, autoregressive models, and variational autoencoders. We highlight the iterative numerical techniques underpinning these models, emphasizing potential advancements via methods like warm restarts and adaptive solvers. Empirical results substantiate the benefits of this unified perspective in terms of accuracy and convergence speed.

Diffusion models.Diffusion models, especially denoising diffusion probabilistic models (DDPMs), capture a noise process and its reverse (denoising) trajectory. While score matching with Langevin dynamics models (SMLDs) is relevant, our focus is primarily on DDPMs for their simpler setup. These models transition from complex pixel-space distributions to more tractable Gaussian distributions. Notably, increasing iterations can enhance the generative quality of DDPMs, a connection we wish to deepen. This procedure can be seen as instantiating an iteration procedure, where iterations are modified as in the methods found in [20]. This operator setting and iterative interpretation is described in detail in Appendix B.1.

To empirically explore convergence with iterations in diffusion models, we train 10 different DDPMs with 100-1000 iterations and analyze their training dynamics and perceptual quality. Figure 2 reveals that increasing timesteps improves FID scores of generated images. Additionally, Figure 2 demonstrates a consistent decrease in both training and test loss with more time steps, attributed to the diminished area under the expected KL divergence curve over time (Figure 8). Notably, FID scores decline beyond the point of test loss convergence, stabilizing after approximately 150,000 steps (Figure 8). This behavior indicates robust convergence with increasing iterations.

AlphaFold.AlphaFold, a revolutionary protein structure prediction model, takes amino acid sequences and predicts their three-dimensional structure. While the model's intricacies can be found in [17], our primary interest lies in mapping AlphaFold within the operator learning context. Considering an input amino acid sequence, it undergoes processing to yield a multiple sequence alignment (MSA) and a pairwise feature representation. These data are subsequently fed into Evoformers and Structure Modules, iteratively refining the protein's predicted structure. We can think of the output of the Evoformer model as pair of functions lying in some discretized Banach space, while the Structure Modules of AlphaFold can be thought of as being operators over a space of matrices. This is described in detail in Appendix B.2.

To empirically explore the convergence behavior of AlphaFold as a function of iterations, we applied AlphaFold-Multimer across a range of 0-20 recycles on each of the 29 monomers using ground truth targets from CASP15. Figure 3 presents the summarized results, which show that while on average the GDT scores and RMSD improve with AlphaFold-Multimer, not all individual targets consistently converge, as depicted in Figures 4 and 5. Given that AlphaFold lacks a convergence constraint in its training, its predictions can exhibit variability across iterations.
Graph Neural Networks.Graph neural networks (GNNs) excel in managing graph-structured data by harnessing a differentiable message-passing mechanism. This enables the network to assimilate information from neighboring nodes to enhance their representations. We can think of the feature spaces as being Banach spaces of functions, which are discretized according to some grid. The GNN architecture can be thought of as being an operator acting on the direct sum of the Banach spaces, where the underlying geometric structure of the graph determines how the operator combines information through the topological information of the graph. A detailed description is given in Appendix A.3, where theoretical guarantees for the convergence of the iterations are given, and Appendix B.3.

Neural Integral Equations.Neural Integral Equations (NIEs), and their variant Attentional Neural Integral Equations (ANIEs), draw inspiration from integral equations. Here, an integral operator, determined by a neural network, plays a pivotal role.

Denoting the integrand of the integral operator as \(G_{\theta}\) within an NIE, the equation becomes:

\[\mathbf{y}=f(\mathbf{y},\mathbf{x},t)+\int_{\Omega\times[0,1]}G_{\theta}( \mathbf{y},\mathbf{x},\mathbf{z},t,s)d\mathbf{z}ds\]

To solve such integral equations, one very often uses iterative methods, as done in [23] and the training of the NIE model consists in finding the parameters \(\theta\) such that the solutions of the corresponding integral equations model the given data. A more detailed discussion of this model is given in Appendix B.4.

## 5 Experiments

In this section, we showcase experiments highlighting the advantages of explicit iterations. We introduce a new GNN architecture based on Picard iteration and enhance vision transformers with

Figure 2: **Left and Middle**: Losses always decrease with more iterations in DDPMs. Training is stable and overfitting never occurs. EMA smoothing with \(\alpha=0.1\) is used for the loss curves to make the differences clearer. **Right**: DDPMs show FID and loss improves with an increased number of iterations on CIFAR-10. The number of iterations represent the denoising steps during training and inference. All diffusion models, UNets of identical architecture, are trained on CIFAR-10’s training dataset with 64-image batches.
Picard iteration.

### PIGN: Picard Iterative Graph Neural Network

To showcase the benefits of explicit iterations in GNNs, we developed **P**icard **I**teration **G**raph neural **N**etwork (PIGN), a GNN that applies Picard iterations for message passing. We evaluate PIGN against state-of-the-art GNN methods and another iterative approach called IterGNN [THG\({}^{+}\)20] on node classification tasks.

GNNs can suffer from over-smoothing and over-squashing, limiting their ability to capture long-range dependencies in graphs [NHN\({}^{+}\)23]. We assess model performance on noisy citation graphs (Cora and CiteSeer) with added drop-in noise. Drop-in noise involves increasing a percentage \(p\) of the bag-of-words feature values, hindering classification. We also evaluate on a long-range benchmark (LRGB) for graph learning [DRG\({}^{+}\)22].

Table 5 shows PIGN substantially improves accuracy over baselines on noisy citation graphs. The explicit iterative process enhances robustness. Table 1 illustrates PIGN outperforms prior iterative and non-iterative GNNs on the long-range LRGB benchmark, using various standard architectures. Applying Picard iterations enables modeling longer-range interactions.

The PIGN experiments demonstrate the benefits of explicit iterative operator learning. Targeting weaknesses of standard GNN training, PIGN effectively handles noise and long-range dependencies. A theoretical study of convergence guarantees is given in Appendix A.3.

### Enhancing Transformers with Picard Iteration

We hypothesize that many neural network frameworks can benefit from Picard iterations. Here, we empirically explore adding iterations to Vision Transformers. Specifically, we demonstrate the

\begin{table}
\begin{tabular}{l c c c} \hline \hline  & **GCN**[15] & **GAT**[15] & **GraphSAGE**[15] \\ \hline w/o iterations & \(0.1510\pm 0.0029\) & \(0.1204\pm 0.0127\) & \(0.3015\pm 0.0032\) \\ IterGNN & \(0.1736\pm 0.0311\) & \(0.1099\pm 0.0459\) & \(0.1816\pm 0.0014\) \\ PIGN (Ours) & \(\mathbf{0.1831\pm 0.0038}\) & \(\mathbf{0.1706\pm 0.0046}\) & \(\mathbf{0.3560\pm 0.0037}\) \\ \hline \hline \end{tabular}
\end{table}
Table 1: F1 scores of different models on the standard test split of the LRGB PascalVOC-SP dataset. Rows refer to model frameworks and columns are GNN backbone layers. A budget of 500k trainable parameters is set for each model. Each model is run on the same set of 5 random seeds. The mean and standard deviation are reported. For IterGNN with GAT backbone, two of the runs keep producing exploding loss so the reported statistics only include three runs.

Figure 3: On average, additional iterations enhance AlphaFold-Multimer’s performance, though the model doesn not invariably converge with more iterations. Target-specific trends can be seen in Figures 4 and 5.
benefits of explicit Picard iteration in transformer models on the task of solving the Navier-Stokes partial differential equation (PDE) as well as self-supervised masked prediction of images. We evaluate various Vision Transformer (ViT) architectures [1] along with Attentional Neural Integral Equations (ANIE) [13].

For each model, we perform training and evaluation with different numbers of Picard iterations as described in Section 3.2. We empirically observe improved performance with more iterations for all models, since additional steps help better approximate solutions to the operator equations.

Table 2 shows lower mean squared error on the PDE task for Vision Transformers when using up to three iterations compared to the standard single-pass models. Table 3 shows a similar trend for self-supervised masked prediction of images. Finally, Table 4 illustrates that higher numbers of iterations in ANIE solvers consistently reduces error. We observe in our experiments across several transformer-based models and datasets that, generally, more iterations improve performance.

Overall, these experiments highlight the benefits of explicit iterative operator learning. For transformer-based architectures, repeating model application enhances convergence to desired solutions. Our unified perspective enables analyzing and improving networks from across domains. A theoretical study of the convergence guarantees of the iterations is given in Appendix A.2.

## 6 Discussion

We introduced an iterative operator learning framework in neural networks, drawing connections between deep learning and numerical analysis. Viewing networks as operators and employing techniques like Picard iteration, we established convergence guarantees. Our empirical results, exemplified by PIGN--an iterative GNN, as well as an iterative vision transformer, underscore the benefits of explicit iterations in modern architectures.

For future work, a deeper analysis is crucial to pinpoint the conditions necessary for convergence

\begin{table}
\begin{tabular}{l c c c} \hline \hline
**Model** & \(N_{\text{iter}}=1\) & \(N_{\text{iter}}=2\) & \(N_{\text{iter}}=3\) \\ \hline ViT & \(0.2472\pm 0.0026\) & \(0.2121\pm 0.0063\) & \(\mathbf{0.0691\pm 0.0024}\) \\ ViTsmall & \(0.2471\pm 0.0025\) & \(0.1672\pm 0.0087\) & \(\mathbf{0.0648\pm 0.0022}\) \\ ViTparallel & \(0.2474\pm 0.0027\) & \(0.2172\pm 0.0066\) & \(\mathbf{0.2079\pm 0.0194}\) \\ ViT3D & \(0.2512\pm 0.0082\) & \(\mathbf{0.2237\pm 0.0196}\) & \(0.2529\pm 00.0079\) \\ \hline \hline \end{tabular}
\end{table}
Table 2: ViT models used to solve a PDE (Navier-Stokes). The mean squared error is reported for each model as the number of iterations varies. A single iteration indicates the baseline ViT model. Higher iterations perform better than the regular ViT (\(N_{\text{iter}}=1\)).
and stability within our iterative paradigm. There remain unanswered theoretical elements about dynamics and generalization. Designing network architectures inherently tailored for iterative processes might allow for a more effective utilization of insights from numerical analysis. We are also intrigued by the potential of adaptive solvers that modify the operator during training, as these could offer notable advantages in both efficiency and flexibility.

In summation, this work shines a light on the synergies between deep learning and numerical analysis, suggesting that the operator-centric viewpoint could foster future innovations in the theory and practical applications of deep neural networks.

\begin{table}
\begin{tabular}{l c c c c c} \hline \hline
**Model size** & \(N_{\text{iter}}=1\) & \(N_{\text{iter}}=2\) & \(N_{\text{iter}}=4\) & \(N_{\text{iter}}=6\) & \(N_{\text{iter}}=8\) \\ \hline \(1H|1B\) & \(0.0564\pm 0.0070\) & \(0.0474\pm 0.0065\) & \(0.0448\pm 0.0062\) & \(0.0446\pm 0.0065\) & \(\mathbf{0.0442\pm 0.0065}\) \\ \(4H|1B\) & \(0.0610\pm 0.0078\) & \(0.0516\pm 0.0083\) & \(0.0512\pm 0.0070\) & \(0.0480\pm 0.0066\) & \(\mathbf{0.0478\pm 0.0066}\) \\ \(2H|2B\) & \(0.0476\pm 0.0065\) & \(0.0465\pm 0.0067\) & \(0.0458\pm 0.0067\) & \(0.0451\pm 0.0064\) & \(\mathbf{0.0439\pm 0.0062}\) \\ \(4H|4B\) & \(0.0458\pm 0.0062\) & \(0.0461\pm 0.0065\) & \(0.0453\pm 0.0063\) & \(0.0453\pm 0.0061\) & \(\mathbf{0.0445\pm 0.0059}\) \\ \hline \hline \end{tabular}
\end{table}
Table 4: Performance of ANIE on a PDE (Navier-Stokes) as the number of iterations of the integral equation solver varies, and for different sizes of architecture. Here \(H\) indicates the number of heads and \(B\) indicates the number of blocks (layers). A single iteration means that the integral operator is applied once. As the number of iterations of the solver increases, the performance of the model in terms of mean squared error improves.

\begin{table}
\begin{tabular}{l c c c} \hline \hline
**Model** & \(N_{\text{iter}}=1\) & \(N_{\text{iter}}=2\) & \(N_{\text{iter}}=3\) \\ \hline ViT (MSE) & \(0.0126\pm 0.0006\) & \(\mathbf{0.0121\pm 0.0006}\) & \(0.0122\pm 0.0006\) \\ ViT (FID) & \(20.0433\) & \(20.0212\) & \(\mathbf{19.2956}\) \\ \hline \hline \end{tabular}
\end{table}
Table 3: ViT models trained with a pixel dropout reconstruction objective on CIFAR-10. The ViT architecture contains 12 encoder layers, 4 decoder layers, 3 attention heads in both the encoder and decoder. The embedding dimension and patch size are 192 and 2. The employed loss is \(\text{MSE}((1-\lambda)T(x_{i})+\lambda x_{i},y)\), computed on the final iteration \(N_{\text{iter}}=i\). Images are altered by blacking out 75% of pixels. During inference, iterative solutions are defined as \(x_{i+1}=(1-\lambda)T(x_{i})+\lambda x_{i}\), for \(i\in\{0,1,\ldots N\}\). Here, \(N=2\) and \(\lambda=1/2\).