Gridded Transformer Neural Processes for Large Unstructured Spatio-Temporal Data

## Abstract

Abstract Many important problems require modelling large-scale spatio-temporal datasets, with one prevalent example being weather forecasting.
Recently, transformer-based approaches have shown great promise in a range of weather forecasting problems. However, these have mostly focused on gridded data sources, neglecting the wealth of unstructured, off-the-grid data from observational measurements such as those at weather stations.
A promising family of models suitable for such tasks are neural processes (NPs), notably the family of transformer neural processes (TNPs). Although TNPs have shown promise on small spatio-temporal datasets, they are unable to scale to the quantities of data used by state-of-the-art weather and climate models. This limitation stems from their lack of efficient attention mechanisms.
We address this shortcoming through the introduction of gridded pseudo-token TNPs which employ specialised encoders and decoders to handle unstructured observations and utilise a processor containing gridded pseudo-tokens that leverage efficient attention mechanisms.
Our method consistently outperforms a range of strong baselines on various synthetic and real-world regression tasks involving large-scale data, while maintaining competitive computational efficiency. The real-life experiments are performed on weather data, demonstrating the potential of our approach to bring performance and computational benefits when applied at scale in a weather modelling pipeline.

## 1 Introduction

Many spatio-temporal modelling problems are being transformed by the proliferation of data from *in situ* sensors, remote observations, and scientific computing models. The opportunities presented have led to a surge of interest from the machine learning community to develop new tools and models to support these efforts.
One prominent example of this transformation is in medium-range weather and environmental forecasting where a new generation of machine learning models have improved performance and reduced computational costs, including: Aurora (Bodnar et al., ); GraphCast (Lam et al., ); GenCast (Price et al., ); PanguWeather (Bi et al., ); ClimaX (Nguyen et al., ), FuXi (Chen et al., ); and FengWu (Chen et al., ). All of these models operate on sets of environmental variables which are regularly structured in space and time, allowing them to leverage architectures developed in the vision and language community such as the Vision Transformer (ViT; Dosovitskiy et al. ), Swin Transformer (Liu et al., ), and Perceiver (Jaegle et al., ).
These models are trained and deployed on data arising from computationally intensive scientific simulation and analysis techniques, which integrate observational data (e.g. from weather stations, ships, buoys, radiosondes, etc.) with simulation data to provide the best estimate of the atmosphere’s state. They are not currently trained on observational measurements directly.

We are now on the cusp of a second generation of these models which will also ingest unstructured observational, alongside analysis data, or which will replace the need for analysis data entirely. This second generation of forecasting systems will further improve accuracy and reduce computational costs.
Thus far, there have only been a handful of early attempts to tackle this problem (Vaughan et al., ; McNally et al., ; Xu et al., ; Xiao et al., ) and it is unclear what architectures would best support this setting. We therefore consider it fertile ground for impactful research.

Many problems considered can be formulated as repeatedly performing predictive inference conditioned on an ever-changing, large set of observations.
Posed with the question of what framework is suitable for handling spatio-temporal prediction problems containing both structured and unstructured data, we turn towards neural processes (NPs; Garnelo et al. ; ): a family of meta-learning models that are able to map from datasets of arbitrary size and structure to predictions over outputs at arbitrary input locations. NPs support a probabilistic treatment of observations, making them capable of outputting uncertainty estimates that are crucial in, for example, forecasting systems. Moreover, NPs are a flexible framework—unlike other models that only perform forecasting, they are able to solve more general state estimation problems, including forecasting, data fusion, data interpolation and data assimilation.
Although early versions of NPs were severely limited, a string of recent developments (Kim et al., ; Gordon et al., ; Nguyen & Grover, ; Ashman et al., ; Feng et al., ; Bruinsma et al., ; Ashman et al., ) have significantly improved their effectiveness, particularly for small-scale spatio-temporal regression problems. Notably, the family of transformer NPs (TNPs; Nguyen & Grover ; Feng et al. ; Ashman et al. ), which use transformer-based architectures as the computational backbone, have demonstrated impressive performance on a range of tasks. However, unlike the aforementioned large-scale environmental prediction models, TNPs are yet to fully take advantage of efficient attention mechanisms. As a result of the quadratic computational complexity associated with full attention, they have been unable to scale to complex spatio-temporal datasets. This is because such techniques require structured—more specifically, gridded—datasets, and are thus not immediately applicable.

We pursue a straightforward solution: to encode the dataset onto a structured grid before passing it through a transformer-based architecture. We refer to this family of models as *gridded TNPs*, and we outline our core contributions as follows:

- 1.
We develop a novel attention-based grid encoder, based on the ideas of ‘pseudo-tokens’ (Jaegle et al., ; Feng et al., ; Lee et al., ), which we show improves upon the performance of traditional kernel-based interpolation methods.
- 2.
Equipped with the pseudo-token grid encoder, we introduce the use of efficient attention mechanisms into the transformer backbone of the TNP, utilising advancements such as the ViT (Dosovitskiy et al., ) and Swin Transformer (Liu et al., ).
- 3.
We develop an efficient $k$-nearest-neighbour attention-based grid decoder, facilitating the evaluation of predictive distributions at arbitrary spatio-temporal locations. Remarkably, we find that this improves performance over full attention.
- 4.
We empirically evaluate our model on a range of synthetic and real-world spatio-temporal regression tasks, demonstrating both the ability to 1) maintain strong performance on large spatio-temporal datasets and 2) handle multiple sources of unstructured data effectively, all while maintaining a low computational complexity.

## 2 Background

We consider the supervised learning setting, where $\mathcal{X}$, $\mathcal{Y}$ denote the input and output spaces, and $(\mathbf{x},\mathbf{y})\in\mathcal{X}\times\mathcal{Y}$ denotes an input-output pair. We restrict our attention to $\mathcal{X}=\mathbb{R}^{D_{x}}$ and $\mathcal{Y}=\mathbb{R}^{D_{y}}$.
Let $\mathcal{S}=\bigcup_{N=0}^{\infty}(\mathcal{X}\times\mathcal{Y})^{N}$ be a collection of all finite datasets, which includes the empty set $\varnothing$.
We denote a context and target set with $\mathcal{D}_{c},\ \mathcal{D}_{t}\in\mathcal{S}$, where $|\mathcal{D}_{c}|=N_{c}$, $|\mathcal{D}_{t}|=N_{t}$.
Let $\mathbf{X}_{c}\in(\mathcal{X})^{N_{c}}$, $\mathbf{Y}_{c}\in(\mathcal{Y})^{N_{c}}$ be the inputs and corresponding outputs of $\mathcal{D}_{c}$, and let $\mathbf{X}_{t}\in(\mathcal{X})^{N_{t}}$, $\mathbf{Y}_{t}\in(\mathcal{Y})^{N_{t}}$ be defined analogously. We denote a single task as $\xi=(\mathcal{D}_{c},\mathcal{D}_{t})=((\mathbf{X}_{c},\mathbf{Y}_{c}),(
\mathbf{X}_{t},\mathbf{Y}_{t}))$.

### 2.1 Neural Processes

Neural processes (NPs; Garnelo et al. ; ) can be viewed as neural-network-based mappings from context sets $\mathcal{D}_{c}$ to predictive distributions at target locations $\mathbf{X}_{t}$, $p(\,\cdot\mid\mathbf{X}_{t},\mathcal{D}_{c})$.
In this work, we restrict our attention to conditional NPs (CNPs; Garnelo et al. ), which only target marginal predictive distributions by assuming that the predictive densities factorise: $p(\mathbf{Y}_{t}|\mathbf{X}_{t},\mathcal{D}_{c})=\prod_{n=1}^{N_{t}}p(\mathbf{
y}_{t,n}|\mathbf{x}_{t,n},\mathcal{D}_{c})$. We denote all parameters of a CNP by $\theta$.
CNPs are trained in a meta-learning fashion, in which the expected predictive log-probability is maximised:

$$ $\textstyle\theta_{\text{ML}}=\operatorname*{arg\,max}_{\theta}\mathcal{L}_{ \text{ML}}(\theta)\quad\text{where}\quad\mathcal{L}_{\text{ML}}(\theta)= \mathbb{E}_{p(\xi)}\big{[}\sum_{n=1}^{N_{t}}\log p_{\theta}(\mathbf{y}_{t,n}| \mathbf{x}_{t,n},\mathcal{D}_{c})\big{]}.$ (1) $$

For real-world datasets, we only have access to a finite number of tasks for training and so approximate this expectation with an average over tasks. The global maximum is achieved if and only if the model recovers the ground-truth predictive distributions (Proposition 3.26 by Bruinsma, ).

#### 2.1.1 A Unifying Construction

The many variants of CNPs differ in their construction of the predictive distribution $p(\cdot|\mathbf{x}_{t},\mathcal{D}_{c})$. Here, we introduce a construction of CNPs that generalises all variants. CNPs are formed from three components: the *encoder*, the *processor*, and the *decoder*. The encoder, $e\colon\mathcal{X}\times\mathcal{Y}\rightarrow\mathcal{Z}$, first encodes each $(\mathbf{x}_{c,n},\mathbf{y}_{c,n})\in\mathcal{D}_{c}$ into some latent representation, or *token*, $\mathbf{z}_{c,n}\in\mathcal{Z}$. The processor, $\rho\colon\big{(}\bigcup_{n=0}^{\infty}\mathcal{Z}^{n}\big{)}\times\mathcal{X}
\rightarrow\mathcal{Z}$, processes the set of context tokens $e(\mathcal{D}_{c})=\{e(\mathbf{x}_{c,n},\mathbf{y}_{c,n})\}_{n}$ together with the target input $\mathbf{x}_{t}$ to obtain a target dependent token, $\mathbf{z}_{t}\in\mathcal{Z}$.(^1^11The space of target dependent tokens does not need to be the same as that of context tokens—we have used $\mathcal{Z}$ in both cases for simplicity. It is also possible for $\mathcal{Z}$ to be the product of multiple spaces, e.g. $\mathcal{Z}=\mathcal{Z}_{token}\times\mathcal{X}$ where we retain information about the input locations.) Finally, the decoder, $d\colon\mathcal{Z}\rightarrow\mathcal{P}_{\mathcal{Y}}$, maps from the target token to the predictive distribution over the output at that target location. Here, $\mathcal{P}_{\mathcal{Y}}$ denotes the set of distributions over $\mathcal{Y}$. When constructing the processor, it is important to ensure that
(a) it operates on the context set of tokens in a permutation invariant manner, and
(b) it is able to operate on multiple target inputs in a single forward pass, which allows for computationally efficient training and inference.
We highlight the order of these three components in Figure 1, and demonstrate how a selection of CNPs can be constructed according to this in Appendix A.

Figure: Figure 1: A unifying construction of CNPs, with $\mathcal{D}_{c}=\{(\mathbf{x}_{c,n},\mathbf{y}_{c,n})\}_{n}$ and $\mathbf{z}_{c,n}=e(\mathbf{x}_{c,n},\mathbf{y}_{c,n})$.

### 2.2 Transformers and Transformer Neural Processes

A useful perspective of transformers is that of set functions (Lee et al., ), making them suitable in the construction of the processor of a CNP, resulting in the family of TNPs. In this section, we provide a brief overview of transformers and their use in TNPs.

#### 2.2.1 Self-Attention and Cross-Attention

Broadly speaking, transformer-based architectures consist of two operations: multi-head self-attention (MHSA) and multi-head cross-attention (MHCA). Informally, the MHSA operation can be understood as updating a set of tokens using the same set, whereas the MHCA operation can be understood as updating one set of tokens using a different set. More formally, let $\mathbf{Z}\in\mathbb{R}^{N\times D_{z}}$ denote a set of $N$ $D_{z}$-dimensional input tokens. The MHSA operation updates this set of tokens as

$$ $\textstyle\mathbf{z}_{n}\leftarrow\operatorname{cat}\Big{(}\Big{\{}\sum_{m=1}^ {N}\alpha_{h}(\mathbf{z}_{n},\mathbf{z}_{m}){\mathbf{z}_{m}}^{T}\mathbf{W}_{V, h}\Big{\}}_{h=1}^{H}\Big{)}\mathbf{W}_{O}\quad\forall\ n=1,\ldots,N.$ (2) $$

Here, $\mathbf{W}_{V,h}\in\mathbb{R}^{D_{z}\times D_{V}}$ and $\mathbf{W}_{O}\in\mathbb{R}^{HD_{V}\times D_{z}}$ are the value and projection weight matrices, where $H$ denotes the number of ‘heads’, and $\alpha_{h}$ is the attention mechanism. This is most often a softmax-normalised transformed inner-product between pairs of tokens: $\alpha_{h}(\mathbf{z}_{n},\mathbf{z}_{m})=\operatorname{softmax}(\{\mathbf{z}_
{n}^{T}\mathbf{W}_{Q,h}\mathbf{W}_{K,h}^{T}\mathbf{z}_{m}\}_{m=1}^{N})_{m}$, where $\mathbf{W}_{Q,h}\in\mathbb{R}^{D_{z}\times D_{QK}}$ and $\mathbf{W}_{K,h}\in\mathbb{R}^{D_{z}\times D_{QK}}$ are the query and key matrices. The MHCA operation updates one set of tokens, $\mathbf{Z}_{1}\in\mathbb{R}^{N_{1}\times D_{z}}$, using another set of tokens, $\mathbf{Z}_{2}\in\mathbb{R}^{N_{2}\times D_{z}}$, in a similar manner:

$$ $\textstyle\mathbf{z}_{1,n}\leftarrow\operatorname{cat}\Big{(}\Big{\{}\sum_{m=1 }^{N_{2}}\alpha_{h}(\mathbf{z}_{1,n},\mathbf{z}_{2,m}){\mathbf{z}_{2,m}}^{T} \mathbf{W}_{V,h}\Big{\}}_{h=1}^{H}\Big{)}\mathbf{W}_{O}\quad\forall\ n=1, \ldots,N_{1}.$ (3) $$

MHSA and MHCA operations are used in combination with layer-normalisation operations and point-wise MLPs to obtain MHSA and MHCA blocks. Unless stated otherwise, we shall adopt the order used by Vaswani et al. () which we detail in Appendix E.

#### 2.2.2 Pseudo-Token-Based Transformers

Pseudo-token-based transformers, first introduced by Jaegle et al. () with the Perceiver, remedy the quadratic computational complexity induced by the standard transformer by condensing the set of $N$ tokens, $\mathbf{Z}\in\mathbb{R}^{N\times D_{z}}$, into a smaller set of $M\ll N$ ‘pseudo-tokens’, $\mathbf{U}\in\mathbb{R}^{M\times D_{z}}$, using MHCA operations.(^2^22There are strong similarities between the use of pseudo-tokens in transformers and the use of inducing points in sparse Gaussian processes: both can be interpreted as condensing the original dataset into a smaller, ‘summary dataset’.) These pseudo-tokens are then processed instead of operating on the original set directly, reducing the computational complexity from $\mathcal{O}(N^{2})$—the cost of performing MHSA operations on the original set—to $\mathcal{O}(NM+M^{2})$—the cost of performing MHCA operations between the original set and pseudo-tokens, followed by MHSA operations on the pseudo-tokens.

#### 2.2.3 Transformer Neural Processes

Transformer neural processes (TNPs) use transformer-based architectures as the processor in the CNP construction described in Section 2.1.1. First, each context point $(\mathbf{x}_{c,n},\mathbf{y}_{c,n})\in\mathcal{D}_{c}$ and target input $\mathbf{x}_{t,n}\in\mathbf{X}_{t}$ is encoded to obtain an initial set of context and target tokens, $\mathbf{Z}^{0}_{c}\in\mathbb{R}^{N_{c}\times D_{z}}$ and $\mathbf{Z}^{0}_{t}\in\mathbb{R}^{N_{t}\times D_{z}}$.
Transformers are then used to process the union $\mathbf{Z}^{0}=\mathbf{Z}^{0}_{c}\cup\mathbf{Z}^{0}_{t}$ using a series of MHSA and MHCA operations, keeping only the output tokens corresponding to the target inputs. These processed target tokens are then mapped to predictive distributions using the decoder. The specific transformer-based architecture is unique to each TNP variant, with each generally consisting of MHSA operations acting on the context tokens—or pseudo-token representation of the context—and MHCA operations acting to update the target tokens given the context tokens. We provide an illustrative diagram of two popular TNP variants in Appendix A: the regular TNP (Nguyen & Grover, ) and
the induced set transformer NP (ISTNP; Lee et al. ).

The application of TNPs to large datasets is impeded by the use of MHSA and MHCA operations acting on the entire set of context and target tokens. Even when pseudo-tokens are used in pseudo-token-based TNPs (PT-TNPs), the number of pseudo-tokens $M$ required for accurate predictive inference generally scales with the complexity of the dataset, and the MHCA operations between the pseudo-tokens and the context and target tokens quickly becomes prohibitive as the size of the dataset increases. This motivates the use of efficient attention mechanisms, as used in the ViT (Dosovitskiy et al., ) and Swin Transformer (Liu et al., ).

## 3 Related Work

### 3.1 Transformers for Point Cloud Data

A closely related research area is point cloud data modelling (Tychola et al., ), with transformer-based architectures employed for a variety of tasks (Lu et al., ). Given the large amount of data used in some works, the use of efficient attention mechanisms has been considered. A number of notable approaches use voxelisation, whereby unstructured point clouds are encoded onto a structured grid (Mao et al., ; Zhang et al., ), prior to employing efficient architectures that operate on grids. Both aforementioned works voxelise the unstructured inputs through rasterisation, but implement efficient self-attention mechanisms such as Swin Transformer.
Our pseudo-token grid encoder is heavily inspired by the voxel-based set attention (VSA) introduced by He et al. (). Similar to our approach, VSA cross-attends local neighbourhoods of unstructured tokens onto a structured grid of pseudo-tokens. However, VSA uses the same initial pseudo-token values for all grid locations, as the ‘cloud’ in which points exist are not equipped with unobserved information. Their use of a fixed set of initial values manifests a specific implementation in which all tokens attend to the same set of pseudo-tokens.
In contrast, for many spatio-temporal problems there may exist potentially unobserved, fixed topographical information such as elevation, land use and soil type, hence we employ different initial pseudo-token values for each grid location which are capable of capturing this.

### 3.2 Models for Structured Weather Data

The motivation behind our method is to develop models that are able to scale to massive spatio-temporal datasets. In recent years, doing so has garnered a significant amount of interest from the research community, and several models have been developed—some of the most successful of which target weather data. However, the majority of these methods operate on structured, gridded data, and focus solely on forecasting.(^3^33That is, to predict the entire gridded state at time $t+1$ given a history of previous gridded state(s).)
Models such as Aurora (Bodnar et al., ), GraphCast (Lam et al., ), GenCast (Price et al., ), Pangu (Bi et al., ), FuXi (Chen et al., ) and FengWu (Chen et al., ) share a similar encoder-processor-decoder to our construction of CNPs presented in Section 2.1, in which the input grid is projected onto some latent space which is then transformed using an efficient form of information propagation (e.g. sparse attention with Swin Transformer (Bodnar et al., ; Bi et al., ), message passing with GNNs (Price et al., )) before being projected back onto the grid at the output.(^4^44It is worth noting that despite transformers being associated with fully-connected graphs, sparse attention mechanisms such as Swin Transformer can also be interpreted as a form of sparsely connected GNN with staggered updates.)

Although yielding impressive performance, these methods are inherently constrained in their application—forecasting with structured data—with predictions sharing the same locations as the inputs at each time point.
The distinct advantage of NPs is their ability to model stochastic processes, which can be evaluated at any target location, and their ability to flexibly condition on unstructured data. Together, these enable NPs to tackle a strictly larger class of tasks.
Indeed, our work reflects a more general trend towards more flexible methods for spatio-temporal data, particularly those that are able to deal with unstructured data. Perhaps the most relevant works to ours are Aardvark (Vaughan et al., ) and FuXi-DA (Xu et al., ), both end-to-end weather prediction models that handle both unstructured and structured data.
Aardvark employs kernel interpolation, also known as a SetConv (Gordon et al., ), to move unstructured data onto a grid, followed by a ViT which processes this grid and outputs gridded predictions, whereas FuXi-DA simply averages observations within each grid cell. We provide an extensive comparison between the kernel interpolation approach to structuring data and our pseudo-token grid encoder in Section 5, demonstrating significantly better performance.
Another relevant example is Lessig et al. (), a task-agnostic stochastic model of atmospheric dynamics, trained to predict randomly masked or distorted tokens. While the task agnosticism of this approach shares similarities with that of NPs, it is unclear how this approach extends to settings in which data can exist at arbitrary spatio-temporal locations, and they are more limited in their approach to model multiple sources of the input data.

## 4 Gridded Transformer Neural Processes

While TNPs have demonstrated promising performance on small to medium-sized datasets, they are unable to scale to the size characteristic of large spatio-temporal datasets. To address this shortcoming, we consider the use of efficient attention mechanisms that have proved effective for gridded spatio-temporal data (Bodnar et al., ; Lam et al., ; Price et al., ; Nguyen et al., ; Bi et al., ). Such methods are not immediately applicable without first structuring the unstructured data. We achieve this by drawing upon methods developed in the point cloud modelling literature—notably the voxel-based set attention (VSA) (He et al., )—and develop the pseudo-token grid encoder: an effective attention-based method for encoding unstructured data onto a grid.

Figure: Figure 2: An illustrative demonstration of the complete gridded TNP pipeline. Following the CNP constrution in Section 2.1.1, we highlight the encoder (blue), processor (red) and decoder (green).
Refer to caption: extracted/5915947/figures/new-pipeline.png

We provide an illustrative diagram of our proposed approach in Figure 2, which decomposes the processor $\rho\colon\big{(}\bigcup_{n=0}^{\infty}\mathcal{Z}^{n}\big{)}\times\mathcal{X}
\rightarrow\mathcal{Z}$ into three parts:
1. the grid encoder, $\rho_{ge}\colon\bigcup_{n=0}^{\infty}\mathcal{Z}^{n}\rightarrow\mathcal{Z}^{M}$, which embeds the set $\{(\mathbf{x}_{c,n},\mathbf{z}_{c,n}\}_{n=1}^{N_{c}}$ into tokens $\{\mathbf{u}_{m}\}_{m=1}^{M}$ at gridded locations $\{\mathbf{v}_{m}\}_{m=1}^{M}$;
2. the grid processor, $\rho_{gp}\colon\mathcal{Z}^{M}\rightarrow\mathcal{Z}^{M}$, which transforms the token values; and
3. the grid decoder, $\rho_{gd}\colon\mathcal{Z}^{M}\times\mathcal{X}\rightarrow\mathcal{Z}$, which maps the tokens onto the location $\mathbf{x}_{t}$.
Here, the latent space $\mathcal{Z}=\mathcal{Z}_{\text{token}}\times\mathcal{X}$. We discuss our choices for each of these components in the remainder of this section.

### 4.1 Grid Encoder: the Pseudo-Token Grid Encoder

Let $\mathbf{z}_{n}\in\mathbb{R}^{D_{z}}$ denote the token representation of input-output pair $(\mathbf{x}_{n},\mathbf{y}_{n})$ after point-wise embedding. We introduce a set of $\prod_{d=1}^{D_{x}}M_{d}$ pseudo-tokens $\mathbf{U}^{0}\in\mathbb{R}^{M_{1}\times\cdots\times M_{D_{x}}\times D_{z}}$ at corresponding locations on the grid $\mathbf{V}\in\mathbb{R}^{M_{1}\times\cdots M_{D_{x}}\times D_{x}}$. That is, pseudo-token $\mathbf{u}_{m_{1},\ldots,m_{D_{x}}}\in\mathbb{R}^{D_{z}}$ is associated with location $\mathbf{v}_{m_{1},\ldots,m_{D_{x}}}\in\mathbb{R}^{D_{x}}$. For ease of reading, we shall replace the product $\prod_{d=1}^{D_{x}}M_{d}$ with $M$ and the indexing notation $m_{1},\ldots,m_{D_{x}}$ with $m$. The pseudo-token grid encoder obtains a pseudo-token representation $\mathbf{U}\in\mathbb{R}^{M\times D_{z}}$ of $\mathcal{D}_{c}$ on the grid $\mathbf{V}$ by cross-attending from each set of tokens $\{\mathbf{z}_{c,n}\}_{n\in\mathfrak{N}(\mathbf{v}_{m};k)}$ to each initial pseudo-token $\mathbf{u}^{0}_{m}$:

$$ $\textstyle\mathbf{u}_{m}\leftarrow\operatorname{MHCA}\left(\mathbf{u}^{0}_{m}, \ \{\mathbf{z}_{c,n}\}_{n\in\mathfrak{N}(\mathbf{v}_{m};k)}\right)\quad\forall m \in M.$ (4) $$

Here, $\mathfrak{N}(\mathbf{v}_{m};k)$ denotes the index set of input locations for which $\mathbf{v}_{m}$ is amongst the $k$ nearest grid locations. In practice, we found that $k=1$ suffices. While this operation seems computationally intensive, there are two tricks we can apply to make it computationally efficient. First, note that provided the pseudo-token grid is regularly spaced, the set of nearest neighbours for all grid locations can be found in $\mathcal{O}(N_{c})$. Second, the operation described in Equation 4 can be performed in parallel for all $m\in M$ by ‘padding’ each set $\{\mathbf{z}_{c,n}\}_{n\in\mathfrak{N}(\mathbf{v}_{m};k)}$ with $\max_{i}|\mathfrak{N}(\mathbf{v}_{i};k)|-|\mathfrak{N}(\mathbf{v}_{m};k)|$ ‘dummy tokens’ and applying appropriate masking, resulting in a computational complexity of $\mathcal{O}\left(M\max_{i}|\mathfrak{N}(\mathbf{v}_{i};k)|\right)$.
Restrictions can be placed on the size of each neighbourhood set to reduce this further. We illustrate the pseudo-token grid encoding in Figure 3.

Figure: Figure 3: An illustrative demonstration of the pseudo-token grid encoder in the 2-D case. To achieve an efficient implementation of cross-attention to the pseudo-token grid, we pad sets of neighbourhood tokens with ‘dummy’ tokens, so that each neighbourhood has the same cardinality.
Refer to caption: extracted/5915947/figures/grid-encoder-new.png

In our experiments, we compare the performance of the pseudo-token grid encoder with simple kernel-interpolation onto a grid, a popular method used by the ConvCNP. We provide a detailed description of this approach in Appendix B.

### 4.2 Grid Processor: Efficient Attention Mechanism-Based Transformers

While we are free to choose any architecture for processing the pseudo-token grid, including CNNs, we focus on transformer-based architectures employing efficient attention mechanisms. Specifically, we consider the use of ViT (Dosovitskiy et al., ) and Swin Transformer (Liu et al., ). Typically, ViT employs *patch encoding* to coarsen a grid of tokens to a coarser grid of tokens, upon which regular MHSA operations are applied. We consider both patch encoding and encoding directly to a smaller grid using the grid encoder, which we collectively refer to when used in gridded TNPs as ViTNPs. The use of Swin Transformer has the advantage of being able to operate with a finer grid of pseudo-tokens, as only local neighbourhoods of tokens attend to each other at any one time. We refer to the use of Swin Transformer in gridded TNPs as Swin-TNP.

### 4.3 Grid Decoder: the Cross-Attention Grid Decoder

In the PT-TNP, *all* pseudo-tokens $\mathbf{U}$ cross-attend to *all* target tokens $\mathbf{Z}_{t}$. This has a computational complexity of $\mathcal{O}(MN_{t})$ which is prohibitive for large $N_{t}$ and $M$. We propose nearest-neighbour cross-attention, in which the set of pseudo-tokens $\{\mathbf{u}_{m}\}_{m\in\tilde{\mathfrak{N}}(\mathbf{x}_{t,n};k)}$ attend to the target token $\mathbf{z}_{t,n}$:

$$ $\textstyle\mathbf{z}_{t,n}\leftarrow\operatorname{MHCA}\big{(}\mathbf{z}^{0}_{ t,n},\ \{\mathbf{u}_{m}\}_{m\in\tilde{\mathfrak{N}}(\mathbf{x}_{t,n};k)}\big{)}.$ (5) $$

Here, $\tilde{\mathfrak{N}}(\mathbf{x}_{t,n};k)$ denotes the index set of the $k$ grid locations that are closest to $\mathbf{x}_{t,n}$.(^5^55More specifically, if the dimension of the target inputs is $d=\text{dim}(x_{t,n})$ for $\forall t,n$, we choose the neighbours within the hypercube defined by the nearest $k^{1/d}$ neighbours in each dimension. When $k^{1/d}$ is not an integer, we take the smallest integer bigger than it (i.e. $\text{ceil}(k^{1/d})$). For more details, please refer to Appendix C.) Assuming equally spaced grid locations, this set can be computed for all target tokens in $\mathcal{O}(N_{t})$, and the operation in Equation 5 can be performed in parallel with a computational complexity of $\mathcal{O}(kN_{t})$. For $k\ll M$, this represents a significant improvement in computational complexity. Further, we found that the use of nearest-neighbour cross-attention improved upon full cross-attention in practice, which we attribute to the inductive biases it introduces being useful for spatio-temporal data.

### 4.4 Handling Multiple Data Modalities

Spatio-temporal datasets often have multiple sources of data, and multi-modal data fusion is an active area of research. It is possible to extend our method to such scenarios when the sources share a common input domain $\mathcal{X}$, the most straightforward of which is to use a different encoder for each modality and apply the pseudo-token grid encoder to the union of tokens. Formally, let $\mathcal{D}_{c}=\bigcup_{s=1}^{S}\mathcal{D}_{c,s}$ denote the context dataset partitioned into $S$ smaller datasets, one for each source. For each source, we obtain $\mathbf{z}_{c,n,s}=e_{s}(\mathbf{x}_{c,n,s},\mathbf{x}_{c,n,s})\ \forall\ (
\mathbf{x}_{c,n,s},\mathbf{y}_{c,n,s})\in\mathcal{D}_{s},\ s\in S$, where $e_{s}\colon\mathcal{X}\times\mathcal{Y}_{s}\rightarrow\mathcal{Z}$ denotes the source-specific point-wise encoder and $\mathcal{Y}_{s}$ denotes the output space for source $s$. We obtain the set of context tokens as $\mathbf{Z}_{c}=\{\{\mathbf{z}_{c,n,s}\}_{n=1}^{N_{c,s}}\}_{s=1}^{S}$. We refer to this approach as the *single* pseudo-token grid encoder.
We also consider a second, more involved approach in which a pseudo-token grid is formed independently for each data modality, with these pseudo-token grids then being merged into a single pseudo-token grid. That is, for each set of source-specific context tokens $\mathbf{Z}_{c,s}=\{e_{s}(\mathbf{x}_{c,n,s},\mathbf{y}_{c,n,s})\}_{n=1}^{N_{c,
s}}$, we obtain a gridded pseudo-token representation $\mathbf{U}_{s}\in\mathbb{R}^{M\times D_{z}}$. A unified gridded pseudo-token representation is obtained by passing their concatenation point-wise through some function. We refer to this approach as the *multi* pseudo-token grid encoder. Note that both methods can be applied analogously to the kernel-interpolation grid encoder.

## 5 Experiments

In this section, we evaluate the performance of our gridded TNPs on synthetic and real-world regression tasks involving datasets with large numbers of datapoints. We demonstrate that gridded TNPs consistently outperform baseline methods, particularly when the difficulty—in terms of both dataset size and complexity—of the predictive inference task increases, while maintaining a low computational complexity. Throughout, we compare the performance of gridded TNPs with two different grid encoders (GEs)—the kernel-interpolation grid encoder (KI-GE) and our pseudo-token grid encoder (PT-GE)—and two different grid processors—Swin Transformer and ViT. We also make comparisons to the following baselines: the CNP (Garnelo et al., ), the PT-TNP using the induced set-transformer architecture (Lee et al., ; Ashman et al., ), and the ConvCNP (Gordon et al., ). We do not compare to the regular ANP or TNP (Kim et al., ; Nguyen & Grover, ) as they are unable to scale to the context set sizes we consider. We provide complete experimental details, including model architectures, datasets and hardware, in Appendices D and E.

### 5.1 Meta-Learning Gaussian Process Regression

We begin with two synthetic 2-D regression tasks with datasets drawn from a Gaussian process (GP) with a squared-exponential (SE) kernel approximated using structure kernel interpolation (SKI) (Wilson & Nickisch, ). Each dataset is a sample of $1.1\times 10^{4}$ datapoints, with input values sampled uniformly from $\mathcal{U}_{[-6,6]}$. For each dataset, we set $N_{c}=1\times 10^{4}$ context and $N_{t}=1\times 10^{3}$ target points. We consider two different setups: one for which the SE kernel lengthscale is fixed to $\ell=0.1$, implying roughly $1.44\times 10^{4}$ ‘wiggles’ in each sampled dataset, and one for which the lengthscale is fixed to $\ell=0.5$, implying roughly $576$ ‘wiggles’. While both setups are non-trivial by the standard of typical NP synthetic experiments, we anticipate the former to be intractable for almost all NP variants—it requires the model to be able to both assimilate finely grained information and modulate its predictions accordingly. In Figure 4, we plot the test log-likelihoods against the time taken to complete a single forward pass (FPT) for a number of gridded TNPs using different grid sizes. We compare the performance to four strong baselines: the PT-TNP with $M=128$ and $M=256$ pseudo-tokens, and the ConvCNP with the same grid sizes as the Swin-TNP.(^6^66Training and inference for all models is performed on a single NVIDIA GeForce RTX 2080 Ti.) A table of results is provided in Section E.1. In Section E.1, we visualise the predictive means of a selection of models for an example dataset with $\ell=0.1$, which demonstrate the superior ability of the gridded TNPs to capture the complexity of the ground-truth dataset.

Figure: Figure 4: Plots comparing the test log-likelihood vs. forward pass time (FPT) for the two synthetic GP datasets. For each model, we show the results for a large and small (transparent) version. The baselines have hatched markers. The grid sizes we consider are $64\times 64$ and $32\times 32$, shown as $64$ and $32$. For the ViTNP models, we include results with and without patch encoding, the former indicated by the $\rightarrow$ symbol in-between the pre- and post-patch-encoded grid sizes. We make use of the following acronyms. KI: kernel-interpolation grid encoding. PT: pseudo-token grid encoding.
Refer to caption: x1.png

### 5.2 Combining Weather Station Observations with Structured Reanalysis

In this experiment, we explore the utility of our model in combining unstructured and sparsely sampled observations with structured, on-the-grid data. We use ERA5 reanalysis data from the European Centre for Medium-Range Weather Forecasts [ECMWF; Hersbach et al. ]. We consider two variables: skin temperature (skt) and 2m temperature (t2m).(^7^77Skin temperature corresponds to the temperature of the Earth’s surface, and is a direct output of standard atmospheric models. In contrast, 2m temperature is found by interpolating between the lowest model temperature (10m) and the skin temperature (Owens & Hewson, ).) We construct each context dataset by combining the t2m at a random subset of $9,957$ weather station locations
(proportion sampled from $\mathcal{U}_{[0,0.3]}$) with a coarsened $180\times 360$ grid—corresponding to a grid spacing of $1^{\circ}$—of skt values.
The target dataset contains the t2m values at all $9,957$ weather station locations. We train on hourly data between 2009-2017, validate on 2018 and provide test metrics on 2019. Training and inference for all models is performed using a single NVIDIA A100 80GB with 32 CPU cores.

In Table 1, we evaluate the performance of six gridded TNPs.
For the Swin-TNP variants, we provide models both with and without the nearest-neighbour cross-attention (NN-CA) mechanism in the decoder. For the Swin-TNP with the PT-GE, we also compare to a model using either solely t2m or skt measurements. Moreover, we provide comparisons to the following baselines: the ConvCNP with a $64\times 128$ grid using a U-Net (Ronneberger et al., ) backbone (and full decoder attention), and the PT-TNP with $M=256$ pseudo-tokens.(^8^88This is the maximum number of pseudo-tokens we could use without running out of memory.) Finally, to study how the models scale with model size, we provide results for Swin-TNP (with PT-GE) and ConvCNP with a grid size of $192\times 384$.

**Table 1: Test log-likelihood ($\boldsymbol{\mathbf{\uparrow}}$) and RMSE ($\boldsymbol{\mathbf{\downarrow}}$) for the t2m station prediction experiment. The standard errors of the log-likelihood are all below $0.010$, and of the RMSE below $0.013$. FPT: forward pass time for a batch size of eight in ms. Params: number of model parameters in units of M.**
| Model | GE | Grid size | Log-lik. $\boldsymbol{\mathbf{\uparrow}}$ | RMSE $\boldsymbol{\mathbf{\downarrow}}$ | FPT | Params |
| --- | --- | --- | --- | --- | --- | --- |
| CNP | - | - | $0.636$ | $3.266$ | $32$ | $0.34$ |
| PT-TNP | - | $M=256$ | $1.344$ | $1.659$ | 230 | $1.57$ |
| ConvCNP (no NN-CA) | SetConv | $64\times 128$ | $1.535$ | $1.252$ | 96 | $9.36$ |
| ViTNP | KI-GE | $48\times 96$ | $1.628$ | $1.197$ | $167$ | $1.14$ |
| ViTNP | PT-GE | $48\times 96$ | $1.704$ | $1.118$ | $181$ | $1.83$ |
| ViTNP | KI-GE | $144\times 288\rightarrow 48\times 96$ | $1.734$ | $1.073$ | $171$ | $1.29$ |
| ViTNP | PT-GE | $144\times 288\rightarrow 48\times 96$ | $1.808$ | $1.021$ | $215$ | $6.69$ |
| Swin-TNP | KI-GE | $64\times 128$ | $1.683$ | $1.157$ | $121$ | $1.14$ |
| Swin-TNP | PT-GE | $64\times 128$ | $\boldsymbol{\mathbf{1.819}}$ | $\boldsymbol{\mathbf{1.006}}$ | $127$ | $2.29$ |
| Swin-TNP (no NN-CA) | KI-GE | $64\times 128$ | $1.544$ | $1.436$ | $137$ | $1.14$ |
| Swin-TNP (no NN-CA) | PT-GE | $64\times 128$ | $1.636$ | $1.273$ | $144$ | $2.29$ |
| Swin-TNP (skt) | PT-GE | $64\times 128$ | $1.427$ | $1.330$ | $123$ | $2.29$ |
| Swin-TNP (t2m) | PT-GE | $64\times 128$ | $1.585$ | $1.599$ | $107$ | $2.24$ |
| ConvCNP | SetConv | $192\times 384$ | $1.689$ | $1.166$ | $74$ | $9.36$ |
| Swin-TNP | PT-GE | $192\times 384$ | $\boldsymbol{\mathbf{2.053}}$ | $\boldsymbol{\mathbf{0.873}}$ | $306$ | $10.67$ |

Figure 5 provides a visual comparison between the t2m predictive errors (i.e. difference between predictive mean and ground truth) of three models: the Swin-TNP with the PT-GE, ConvCNP, and PT-TNP. We show the results for the US, and provide results for the world in Figure 16. The most prominent difference between the predictions is at the stations within central US, where both ConvCNP and PT-TNP tend to underestimate the temperature, whereas the Swin-TNP has relatively small predictive errors. We perform an analysis of the predictive uncertainties in Section E.2.1.

Figure: (a) Swin-TNP error.
Refer to caption: x2.png

### 5.3 Combining Multiple Sources of Unstructured Wind Speed Observations

In our final experiment, we evaluate the utility of our gridded TNPs in modelling sparsely sampled eastward ($\verb|u|$) and northward ($\verb|v|$) components of wind at three different pressure levels: 700hPa, 850hPa and 1000hPa (surface level). Each of the six modalities are obtained from the ERA5 reanalysis dataset (Hersbach et al., ). We perform spatio-temporal interpolation over a latitude / longitude range of $[25^{\circ},\ 49^{\circ}]$ / $[-125^{\circ},-66^{\circ}]$, which corresponds to the contiguous US, spanning four hours. The proportion of the $543,744$ observations used as the context dataset is sampled from $\mathcal{U}_{[0.05,0.25]}$. The target set size is fixed at $135,936$ (i.e. $25\%$ of observations).

We evaluate the performance of eight gridded TNPs: the ViTNP with a $4\times 12\times 30$ grid of pseudo-tokens using the multi grid-encoding method presented in Section 4.4, both the KI-GE and PT-GE, and with and without patch encoding from a $4\times 24\times 60$ grid; and the Swin-TNP with a $4\times 24\times 60$ grid of pseudo-tokens using both multi-modal grid-encoding methods presented in Section 4.4, and both the KI-GE and PT-GE. We provide comparisons to the following baselines: the ConvCNP(^9^99For the ConvCNP, each modality is treated as a different input channel as in Vaughan et al. ().) with a $4\times 24\times 60$ grid using a U-Net backbone, and the PT-TNP with $M=64$ pseudo-tokens.(^10^1010We found this to be the maximum number of pseudo-tokens we could use before running out of memory.) We also provide results for the ConvCNP and Swin-TNP with the PT-GE with larger grid sizes. We visualise the predictive errors for the Swin-TNP, ConvCNP and PT-TNP in Section E.3, which demonstrate the Swin-TNP also has the most accurate predictive uncertainties.

**Table 2: Test log-likelihood ($\boldsymbol{\mathbf{\uparrow}}$) and RMSE ($\boldsymbol{\mathbf{\downarrow}}$) for the the multi-modal wind speed dataset. All grids have a grid size of $4$ in the first dimension (time). The standard errors of the log-likelihoods are all below $0.02$, and of the RMSE below $0.005$. FPT: forward pass time for a batch size of eight in ms. Params: number of model parameters in units of M.**
| Model | GE | Grid size | Log-lik. $\boldsymbol{\mathbf{\uparrow}}$ | RMSE $\boldsymbol{\mathbf{\downarrow}}$ | FPT | Params |
| --- | --- | --- | --- | --- | --- | --- |
| CNP | - | - | $-1.593$ | $2.536$ | $33$ | $0.66$ |
| PT-TNP | - | $M=64$ | $3.988$ | $1.185$ | $166$ | $1.79$ |
| ConvCNP | SetConv | $24\times 60$ | $6.143$ | $0.784$ | $210$ | $14.41$ |
| ViTNP | multi KI-GE | $12\times 30$ | $5.371$ | $0.908$ | $349$ | $1.49$ |
| ViTNP | multi PT-GE | $12\times 30$ | $7.754$ | $0.651$ | $374$ | $2.36$ |
| ViTNP | multi KI-GE | $24\times 60\rightarrow 12\times 30$ | $6.906$ | $0.718$ | $372$ | $1.55$ |
| ViTNP | multi PT-GE | $24\times 60\rightarrow 12\times 30$ | $7.288$ | $0.681$ | $392$ | $3.25$ |
| Swin-TNP | multi KI-GE | $24\times 60$ | $7.603$ | $0.651$ | $355$ | $1.49$ |
| Swin-TNP | multi PT-GE | $24\times 60$ | $\boldsymbol{\mathbf{8.603}}$ | $\boldsymbol{\mathbf{0.577}}$ | $375$ | $3.19$ |
| Swin-TNP | single KI-GE | $24\times 60$ | $7.794$ | $0.642$ | $366$ | $1.39$ |
| Swin-TNP | single PT-GE | $24\times 60$ | $8.073$ | $0.614$ | $364$ | $1.67$ |
| ConvCNP | SetConv | $48\times 120$ | $7.841$ | $0.615$ | $64$(^11^1111The FPT for the larger ConvCNP is smaller than for the smaller ConvCNP as we were forced to use $k=27$ nearest-neighbour grid decoding to avoid out of memory issues.) | $14.41$ |
| Swin-TNP | multi PT-GE | $48\times 120$ | $\boldsymbol{\mathbf{9.383}}$ | $\boldsymbol{\mathbf{0.509}}$ | $369$ | $6.51$ |

### 5.4 Discussion of Results

Across all three experiments, we draw the following conclusions:
1. gridded TNPs significantly outperform all baselines while maintaining relatively low computational complexity;
2. the PT-GE outperforms the KI-GE in almost all gridded TNPs, particularly when the grid is coarse;
3. when modelling multiple data sources, the use of source-specific grid encoders outperforms a single, unified grid encoder;
4. Swin-TNP achieves either comparable or better performance than the ViTNP at a smaller computational cost; and
5. in the grid decoder, nearest-neighbour cross-attention both reduces the computational cost and improves predictive performance relative to full cross-attention.
In Appendix E, we provide additional experimental results showing that these conclusions are robust to changes in model size and dataset composition.

## 6 Conclusion

This paper introduces gridded TNPs, an extension to the family of TNPs which facilitates the use of efficient attention-based transformer architectures such as the ViT and Swin Transformer. Gridded TNPs decompose the computational backbone of TNPs into three distinct components:
1. the *grid encoder*, which moves point-wise data representations onto a structured grid of pseudo-tokens;
2. the *grid processor*, which processes this grid using computationally efficient operations; and
3. the *grid decoder*, which evaluates the processed grid at arbitrary input locations.
In achieving this, we develop the pseudo-token grid encoder—a novel approach to moving unstructured spatio-temporal data onto a grid of pseudo-tokens—and the pseudo-token grid decoder—a computationally efficient approach of evaluating a grid of pseudo-tokens at arbitrary input locations. We compare the performance of a number of gridded TNPs against several strong baselines, demonstrating significantly better performance on large-scale synthetic and real-world regression tasks involving context sets containing over $100,000$ datapoints. This work marks an initial step towards building architectures for modelling large amounts of unstructured spatio-temporal observations. We believe that the methods developed in this paper can be used to both improve and broaden the capabilities of existing machine learning models for spatio-temporal data, including those targeting weather and environmental forecasting. We look forward to pursuing this in future work.

## Acknowledgements

CD is supported by the Cambridge Trust Scholarship. AW acknowledges support from a Turing AI fellowship under grant EP/V025279/1 and the Leverhulme Trust via CFI. RET is supported by Google, Amazon, ARM, Improbable, EPSRC grant EP/T005386/1, and the EPSRC Probabilistic AI Hub (ProbAI, EP/Y028783/1). We thank Anna Vaughaun, Stratis Markou and Wessel Bruinsma for their valuable feedback.

## Appendix A Conditional Neural Processes Architectures

We present below several schematics showing the architectures of different members of the CNP family, and detail how they can be constructed following the universal construction outlined in Section 2.1.1.

##### Original CNP

The first architecture is based on Garnelo et al. () and is the least complex of the CNP variants, using a summation as the permutation-invariant aggregation. The diagram is shown in Figure 6. The encoder of a CNP is an MLP, which maps from each concatenated pair $(\mathbf{x}_{c,n},\mathbf{y}_{c,n})\in\mathcal{D}_{c}$ to some representation $\mathbf{z}_{c,n}\in\mathbb{R}^{D_{z}}$. The processor sums together these representations, and combines the aggregated representation with the target input using $\tau$: $\rho(\{\mathbf{z}_{c,n}\}_{n},\mathbf{x}_{t})=\tau(\sum_{n}\mathbf{z}_{c,n},
\mathbf{x}_{t})$. $\tau$ is often just the concatenation operation. Finally, the decoder consists of another MLP which maps from $\mathbf{z}_{t}=\tau(\sum_{n}\mathbf{z}_{c,n},\mathbf{x}_{t})$ to the parameter space of some distribution over the output space (e.g. Gaussian).

Figure: Figure 6: A diagram illustrating the architecture of the plain CNP (Garnelo et al., ). First, the context set $(\mathbf{x}_{c,n},\mathbf{y}_{c,n})$ and the target tokens $\mathbf{x}_{t,n}$ are encoded using point-wise embeddings. These are fed into the CNP processor, which performs a simple permutation-invariant aggregation of the context tokens. These are then concatenated with the target tokens and fed into the decoder, which outputs the parameters of the specified NP distribution based on the target representation (in this case, mean and variance of a Gaussian).
Refer to caption: extracted/5915947/figures/tnp_architectures/CNP_diagram_v2.png

##### ConvCNP

Another member of the CNP family is the ConvCNP (Gordon et al., ), that embeds sets into function space in order to achieve translation equivariance. This can lead to more efficient training in applications where such an inductive bias is appropriate, such as stationary time-series or spatio-temporal regression tasks. We show in Figure 7 a schematic of the architecture. The encoder of a ConvCNP is simply the identity function. The processor then encodes the discrete function represented by input-output pairs $\{(\mathbf{x}_{c,n},\mathbf{y}_{c,n})\}_{n}$ onto a regular grid using the kernel-interpolation grid encoder (KI-GE). It then processes this grid using a CNN, which is afterwards combined with the target location using a kernel-interpolation grid decoder (KI-GD). Letting $\mathbf{U}=\{\mathbf{u}_{m}\}_{m}$ denote the set of $M$ values on gridded locations $\mathbf{V}=\{\mathbf{v}_{m}\}_{m}$, we can decompose the ConvCNP processor as

$$ $\displaystyle\mathbf{u}_{m}$ $\displaystyle\leftarrow\sum_{n}[1,\ \mathbf{y}_{c,n}]^{T}\psi_{ge}(\mathbf{v}_ {m}-\mathbf{x}_{c,n})$ (6) $\displaystyle\mathbf{U}$ $\displaystyle\leftarrow\operatorname{CNN}(\mathbf{U},\mathbf{V})$ (7) $\displaystyle\mathbf{z}_{t}$ $\displaystyle\leftarrow\sum_{m}\mathbf{u}_{m}\psi_{gd}(\mathbf{x}_{t}-\mathbf{ v}_{m}).$ (8) $$

Here, $\psi_{ge},\ \psi_{gd}\colon\mathbb{R}^{D_{x}}\times\mathbb{R}^{D_{x}}
\rightarrow\mathbb{R}$ denote the KI-GE kernel and KI-GD kernel. As with the CNP, the decoder is an MLP mapping to the parameter space of some distribution over the output space.

Figure: Figure 7: A diagram illustrating the architecture of the ConvCNP (Gordon et al., ). First, the context set $(\mathbf{x}_{c,n},\mathbf{y}_{c,n})$ and the target tokens $\mathbf{x}_{t,n}$ are encoded using point-wise embeddings. These are fed into the ConvCNP processor, which uses a KI-GE to project the tokens into function space. These are then evaluated at discrete locations using a pre-specified resolution, followed by multiple layers of a CNN-based architecture acting upon the discretised signal. To decode at arbitrary locations, a KI-GD is used, giving rise to the target token representation. This is fed into the decoder, which outputs the parameters of the specified NP distribution (in this case, mean and variance of a Gaussian).
Refer to caption: extracted/5915947/figures/tnp_architectures/ConvCNP_diagram_v2.png

##### TNPs

There are a number of different architectures used for the different members of the TNP family. We provide below diagrams for two members mentioned in the main paper, namely the TNP of Nguyen & Grover () and the induced set transformer (ISTNP) of Lee et al. (). For the standard TNP, the encoder consists of an MLP mapping from each $(\mathbf{x}_{c,n},\mathbf{y}_{c,n})\in\mathcal{D}_{c}$ to some representation (token) $\mathbf{z}_{c,n}\in\mathbb{R}^{D_{z}}$. The processor begins by embedding the target location in the same space as the context tokens, giving $\mathbf{z}_{t}\in\mathbb{R}^{D_{z}}$. It then iterates between applying MHSA operations on the set of context tokens, and MHCA operations from the set of context tokens into the target token:

$$ $\left.\begin{matrix}\mathbf{Z}_{c}\leftarrow\operatorname{MHSA}(\mathbf{Z}_{c} )\\ \mathbf{z}_{t}\leftarrow\operatorname{MHCA}(\mathbf{z}_{t};\ \mathbf{Z}_{c}) \end{matrix}\right\}\times L.$ (9) $$

Again, the decoder consists of an MLP mapping from $\mathbf{z}_{t}$ to the parameter space of some distribution over outputs.

Figure: Figure 8: A diagram illustrating the architecture of the TNP (Nguyen & Grover, ). First, the context set $(\mathbf{x}_{c,n},\mathbf{y}_{c,n})$ and the target tokens $\mathbf{x}_{t,n}$ are encoded using point-wise embeddings to obtain the context set representation $\mathbf{Z}_{c}$ and target representation $\mathbf{Z}_{t}$. These are fed into the TNP processor, which takes in the union of $[\mathbf{Z}_{c},\mathbf{Z}_{t}]$ and outputs the token corresponding to the target inputs $\mathbf{Z}^{(L)}_{t}$. At each layer of the processor, the context set representation is first updated through an MHSA layer, which is then used to modulate the target set representation through a MHCA layer between the target set representation from the previous layer and the updated context representation. Finally, the decoder outputs the parameters of the specified NP distribution based on the target representation from the final layer (in this case, mean and variance of a Gaussian).
Refer to caption: extracted/5915947/figures/tnp_architectures/TNP_diagram_v2.png

One of the main limitations of TNPs is the cost of the attention mechanism, which scales quadratically with the number of input tokens. Several works (Feng et al., ; Lee et al., ) addressed this shortcoming by incorporating ideas from the Perceiver-style architecture (Jaegle et al., ) into NPs. The strategy is to introduce a set of $M$ ‘pseudo-tokens’ which act as an information bottleneck between the context and target sets. Provided that $M<<N_{c}$, where $N_{c}$ is the number of context points, this leads to a significant reduction in computational complexity. The architecture we consider in this work is called the induced set transformer NP (ISTNP), and differs from the plain TNP in the calculations performed in the processor. At each layer, the pseudo-token representation is first updated through an MHCA operation from the context set to the pseudo-tokens. The updated representation is then used to modulate the context and target sets separately, through separate MHCA operations:

$$ $\left.\begin{matrix}\mathbf{U}\leftarrow\operatorname{MHCA}(\mathbf{U};\ \mathbf{Z}_{c})\\ \mathbf{z}_{t}\leftarrow\operatorname{MHCA}(\mathbf{z}_{t};\ \mathbf{U})\\ \mathbf{Z}_{c}\leftarrow\operatorname{MHCA}(\mathbf{Z}_{c};\ \mathbf{U})\end{ matrix}\right\}\times L.$ (10) $$

Thus, as apparent in Figure 9, the context and target sets never interact directly, but only through the ‘pseudo-tokens’. The computational cost at each layer reduces from $\mathcal{O}(N_{c}^{2}+N_{c}N_{t})$ in the plain TNP, where $N_{c}$ and $N_{t}$ represent the number of context and target points, respectively, to $\mathcal{O}(M(2N_{c}+N_{t}))$. This is a significant reduction provided that $M<<N_{c}$, resulting in an apparent linear dependency on $N_{c}$. However, in practice, $M$ is not independent of $N_{c}$, with more pseudo-tokens needed as the size of context set increases.

Figure: Figure 9: A diagram illustrating the architecture of the ISTNP (Lee et al., ). As opposed to the regular TNP of Nguyen & Grover (), the ISTNP uses a ‘summarised’ representation of the context set through the use of pseudo-tokens. They are first randomly initialised ($\mathbf{U}^{(0)}$). Then, their representation is updated through cross attention with the context set representation from the previous layer (i.e. at layer $l$: $\mathbf{U}^{(l)}=\text{MHCA}^{(l)}(\mathbf{U}^{(l-1)},\mathbf{Z}_{c}^{(l)})$). This updated set of pseudo-tokens $\mathbf{U}^{(l)}$ is then used to modulate both the context set representation at the current layer through cross-attention (i.e. $\mathbf{Z}_{c}^{(l)}=\text{MHCA}(\mathbf{Z}_{c}^{(l-1)},\mathbf{U}^{(l)})$), as well as the target set representation (i.e. $\mathbf{Z}_{t}^{(l)}=\text{MHCA}(\mathbf{Z}_{t}^{(l-1)},\mathbf{U}^{(l)})$). Thus, the context and target set representations do not interact directly, but only through the pseudo-tokens, which act as a bottleneck of information flow between the two in order to decrease the computational demands of the plain TNP.
Refer to caption: extracted/5915947/figures/tnp_architectures/ISTNP_diagram_v2.png

## Appendix B Kernel-Interpolation Grid Encoder

Let $\mathbf{z}_{n}\in\mathbb{R}^{D_{z}}$ denote the token representation of input-output pair $(\mathbf{x}_{n},\mathbf{y}_{n})$ after point-wise embedding. We introduce the set of grid locations $\mathbf{V}\in\mathbb{R}^{M_{1}\times\cdots M_{D_{x}}\times D_{x}}$. For ease of reading, we shall replace the product $\prod_{d=1}^{D_{x}}M_{d}$ with $M$ and the indexing notation $m_{1},\ldots,m_{D_{x}}$ with $m$.

The kernel-interpolation grid encoder obtains a pseudo-token representation $\mathbf{U}\in\mathbb{R}^{M\times D_{z}}$ of $\mathcal{D}_{c}$ on the grid $\mathbf{V}$ by interpolating from *all* tokens $\{\mathbf{z}_{c,n}\}_{n}$ at corresponding locations $\{\mathbf{x}_{c,n}\}_{n}$ to *all* pseudo-token locations:

$$ $\mathbf{u}_{m}\leftarrow\sum_{n}\mathbf{z}_{c,n}\psi(\mathbf{v}_{m},\mathbf{x} _{c,n})\quad\forall m\in M.$ (11) $$

Here, $\psi\colon\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$ is the kernel used for interpolation, which we take to be the squared-exponential (SE) kernel when $\mathcal{X}=\mathbb{R}^{D_{x}}$:

$$ $\psi_{SE}(\mathbf{v}_{m},\mathbf{x}_{c,n})=\exp\left(-\sum_{d=1}^{D_{x}}\frac{ (x_{c,n,d}-v_{m,d})^{2}}{\ell_{d}^{2}}\right)$ (12) $$

where $\ell_{d}$ denotes the ‘lengthscale’ for dimension $d$. Similar to the pseudo-token grid encoder, we can restrict the kernel-interpolation grid encoder to interpolate only from sets of token $\{\mathbf{z}_{c,n}\}_{n\in\mathfrak{N}(\mathbf{v}_{m};k)}$ for which $\mathbf{v}_{m}$ is amongst the $k$ nearest grid locations. The computational complexity of the kernel-interpolation and pseudo-token grid encoders differ only by a scale factor when this approach is used.

## Appendix C Nearest-Neighbour Cross-Attention in the Grid Decoder

In this section we provide more details, alongside schematics, of our nearest-neighbour cross-attention scheme. As mentioned in the main text, in the grid decoder we only allow a subset of the gridded pseudo-tokens to attend to the target token (those in the vicinity of it). Finding the nearest neighbours is not done in the standard k-nearest neighbours fashion, because we use the same number of nearest neighbours along each dimension of the original data (i.e. latitude and longitude for spatial interpolation; latitude, longitude and time for spatio-temporal interpolation). This is to ensure that we do not introduce specific preferences for any dimension. In the case the data lies on a grid with the same spacing in all its dimensions, the procedure becomes equivalent to k-nearest neighbours.

In practice, we specify the total number of nearest-neighbours we want to use $k$. We then compute the number of nearest-neighbours in each dimension by $k_{\text{dim}}=\text{ceil}(k^{\frac{1}{\text{dim}(x)}})$, where $\text{dim}(x)$ represents the dimensionality of the input. For efficient batching purposes, we tend to choose $k=(2n-1)^{\text{dim}(x)}$, where $n\in\mathbb{N}$ (i.e. $9$ for experiments with latitude-longitude grids, $27$ for experiments with latitude-longitude-time grids). We then find the indices in each dimension of these nearest-neighbours by performing an efficient search that leverages the gridded nature of the data, leading to a computational complexity of $\mathcal{O}(kN_{t})$, with $N_{t}$ the number of target points. When the neighbours go off the grid (i.e. for targets very close to the edges of the grid), we only consider the number of viable (i.e. within the bounds of the grid) neighbours.

##### Example in 2D

This procedure is visualised in Figure 10, where we consider both grids with the same spacing along each dimension, as well as grids with different spacings. We cover both the case of a central target point, as well as a target point closer to the edges of the grid.

Figure: (a) Same / Central
Refer to caption: extracted/5915947/figures/nearest_neighbour/Grid1.png

##### Accounting for non-Euclidean geometry

A lot of our experiments are performed on environmental data, distributed across the Earth. For our purposes, we assume the Earth shows cylindrical geometry, whereby there is no such thing as a grid edge along the longitudinal direction (i.e. a longitude of -180^∘ is the same as 180^∘). This is not the same for latitude, where one extreme corresponds to the North Pole, and the other one to the South Pole. Thus, we would like to allow for the grid to ‘roll’ around the longitudinal direction when computing the nearest neighbours. In this case, there should be no edge target points along the longitudinal dimension. This procedure is graphically depicted in Figure 11 and we use it in our experiments on environmental data.

Figure: Figure 11: Example of nearest-neighbours procedure in $2$D on a $5\times 9$ grid for $9$ nearest neighbours. We allow rolling along the horizontal dimension (e.g. longitude), but do not allow rolling along the vertical one (i.e. latitude). Hence, the neighbours extend on the other side of the grid horizontally, but not vertically. This example corresponds to cylindrical symmetry.
Refer to caption: extracted/5915947/figures/nearest_neighbour/Grid5.png

##### Example for 3D data

We also provide examples of the nearest-neighbours procedure on 3D spaces. These dimensions could represent, for example, latitude, longitude and time, or latitude, longitude and height/pressure levels. In Figure 12 we consider a case where we do not allow for rolling along any dimension, while in Figure 13 we allow for rolling along one of the dimensions.

Figure: (a) Central
Refer to caption: extracted/5915947/figures/nearest_neighbour/Grid_3D_1.png

Figure: (a) Central
Refer to caption: extracted/5915947/figures/nearest_neighbour/Grid_3D_rolling1.png

## Appendix D Hardware specifications

For the smaller synthetic GP regression experiment, we perform training and inference for all models on a single NVIDIA GeForce RTX 2080 Ti GPU with 20 CPU cores. For the other two, larger experiments, we perform training and inference for all models on a single NVIDIA A100 80GB GPU with 32 CPU cores.

## Appendix E Experiment Details

##### Common optimiser details

For all experiments and all models, we use the AdamW optimiser (Loshchilov & Hutter, ) with a fixed learning rate of $5\times 10^{-4}$ and apply gradient clipping to gradients with magnitude greater than 0.5.

##### Common likelihood details

For all experiments and all models, we employ a Gaussian likelihood parameterised by a mean and inverse-softplus variance, i.e. the decoder of each model outputs

$$ $\boldsymbol{\mathbf{\mu}}_{t},\ \log(\exp{\boldsymbol{\mathbf{\sigma}}_{t}^{2} }-1)=d(\mathbf{z}_{t}),\quad p(\ \cdot\mid\mathcal{D}_{c},\mathbf{x}_{t})= \mathcal{N}(\cdot;\boldsymbol{\mathbf{\mu}}_{t},\boldsymbol{\mathbf{\sigma}}_{ t}^{2}).$ (13) $$

For the experiment modelling skin and 2m temperature (skt and t2m) with a richer context, we set a minimum noise level of $\boldsymbol{\mathbf{\sigma}}_{\text{min}}^{2}=0.01$ by parameterising

$$ $\boldsymbol{\mathbf{\mu}}_{t},\ \log(\exp{(\boldsymbol{\mathbf{\sigma}}_{t}- \boldsymbol{\mathbf{\sigma}}_{\text{min}})^{2}}-1)=d(\mathbf{z}_{t}),\quad p( \ \cdot\mid\mathcal{D}_{c},\mathbf{x}_{t})=\mathcal{N}(\cdot;\boldsymbol{ \mathbf{\mu}}_{t},\boldsymbol{\mathbf{\sigma}}_{t}^{2}).$ (14) $$

##### CNP details

For the CNPs, we encode each $(\mathbf{x}_{c,n},\mathbf{y}_{c,n})\in\mathcal{D}_{c}$ in $\mathbb{R}^{D_{z}}$ using an MLP with two-hidden layers of dimension $D_{z}$. We obtain a representation for the entire context set by summing these representations together, $\mathbf{z}_{c}=\sum_{n}\mathbf{z}_{c,n}$, which is then concatenated with the target input $\mathbf{x}_{t}$. The concatenation $[\mathbf{z}_{c},\ \mathbf{x}_{t}]$ is decoded using an MLP with two-hidden layers of dimension $D_{z}$. We use $D_{z}=128$ in all experiments.

##### ConvCNP details

For the ConvCNP model, we use a U-Net architecture (Ronneberger et al., ) for the CNN consisting of 11 layers with input size $C$. Between the five downward layers we apply pooling with size two. For the five upward layers, we use $2C$ input channels and $C$ output channels, as the input channels are formed from the output of the previous layer concatenated with the output of the corresponding downward layer. Between the upward layers we apply linear up-sampling to match the grid size of the downward layer. In all experiments, we use $C=128$, a kernel size of five or nine, and a stride of one. We use SE kernels for the SetConv encoder and SetConv decoder with learnable lengthscales for each input dimension. The grid encoding is modified similarly to the pseudo-token grid encoder, whereby we only interpolate from the set of observations for which each grid point is the closest grid point. Unless otherwise specified, we also modify the grid decoding similarly to the pseudo-token grid decoder, whereby we only interpolate from the $k=3^{D_{x}}$ nearest points on a distance-normalised grid to each target location. We resize the output of the SetConv encoder to dimension $C$ using an MLP with two hidden layers of dimension $C$. We resize the output of the SetConv decoder using an MLP with two hidden layers of dimension $C$.

##### Common transformer details

For each MHSA / MHCA operation, we construct a layer consisting of two residual connections, two layer norm operations, one MLP, together with the MHSA / MHCA operation as follows:

$$ $\displaystyle\widetilde{\mathbf{Z}}$ $\displaystyle\leftarrow\mathbf{Z}+\operatorname{MHSA/MHCA}(\operatorname{layer -norm}_{1}(\mathbf{Z}))$ (15) $\displaystyle\mathbf{Z}$ $\displaystyle\leftarrow\widetilde{\mathbf{Z}}+\operatorname{MLP}(\operatorname {layer-norm}_{2}(\widetilde{\mathbf{Z}})).$ $$

All MHSA / MHCA operations use $H=8$ heads, each with $D_{V}=16$ dimensions. We use a $D_{z}=D_{QK}=128$ throughout.

##### PT-TNP details

We use an induced set transformer (IST) architecture for the PT-TNPs, with each layer consisting of the following set of operations:

$$ $\displaystyle\mathbf{U}$ $\displaystyle\leftarrow\operatorname{MHCA-layer}(\mathbf{U};\mathbf{Z}_{c})$ (16) $\displaystyle\mathbf{z}_{t}$ $\displaystyle\leftarrow\operatorname{MHCA-layer}(\mathbf{z}_{t};\mathbf{U})$ $\displaystyle\mathbf{Z}_{c}$ $\displaystyle\leftarrow\operatorname{MHCA-layer}(\mathbf{Z}_{c};\mathbf{U}).$ $$

In all experiments we use five layers, and encoder / decoder MLPs consisting of two hidden layers of dimension $D_{z}$.

##### ViT details

The ViT architecture consists of optional patch encoding, followed by five MHSA layers. The patch encoding is implemented using a single linear layer. In all experiments we use five layers, and encoder / decoder MLPs consisting of two hidden layers of dimension $D_{z}$.

##### Swin Transformer details

Each layer of the Swin Transformer consists of two MHSA layers applied to each window, and a shifting operation between them. Unless otherwise specified, we use a window size of four and shift size of two for all dimensions (except for the time dimension in the final experiment, as the original grid only has four elements in the time dimension). For the second experiment in which the grid covers the entire globe, we allow the Swin attention masks to ‘roll’ over the longitudinal dimension, allowing the pseudo-tokens near $180^{\circ}$ longitude to attend to those near $-180^{\circ}$ longitude. In all experiments, we use five Swin Transformer layers (10 MHSA layers in total), and encoder / decoder MLPs consisting of two hidden layers of dimension $D_{z}$. We found that the use of a hierarchical Swin Transformer—as used in the original Swin Transformer and models such as Aurora (Bodnar et al., )—did not lead to any improvement in performance.

##### Spherical harmonic embeddings

When modelling input data on the sphere (i.e. the final two experiments), the CNP, PT-TNP, and gridded TNP models first encode the latitude / longitude coordinates using spherical harmonic embeddings following Rußwurm et al. () using 10 Legendre polynomials. We found this to improve performance in all cases.

##### Temporal Fourier embeddings

When modelling input data through time (i.e. the final experiment), the CNP, PT-TNP, and gridded TNP models first encode the temporal coordinates using a Fourier embedding. The time value is originally provided in hours since 1st January 1970. Following Bodnar et al. (), we embed this using a Fourier embedding of the following form:

$$ $\operatorname{Emb}(t)=\big{[}\cos\frac{2\pi t}{\lambda_{i}},\ \sin\frac{2\pi t }{\lambda_{i}}\big{]}\ \text{for}\ 0\leq i<L/2.$ (17) $$

where the $\lambda_{i}$ are log-spaced values between the minimum and maximum wavelength. We set $\lambda_{min}=1$ and $\lambda_{max}=8760$, the number of hours in a year. We use $L=10$.

##### Great-circle distance

For methods using the kernel-interpolation grid encoder (i.e. the ConvCNP and some gridded TNPs), we use the great-circle distance, rather than Euclidean distance, as the input into the kernel when modelling input data on the sphere. The haversine formula determines the great-circle distance between two points $\mathbf{x}_{1}=(\varphi_{1},\lambda_{2})$ and $\mathbf{x}_{2}=(\varphi_{2},\lambda_{2})$, where $\lambda$ and $\varphi$ denote the latitude and longitude, and is given by:

$$ $d(\mathbf{x}_{1},\mathbf{x}_{2})=2r\arcsin\left(\sqrt{\frac{1-\cos(\Delta \varphi)+\cos\varphi_{1}\cdot\cos\varphi_{2}\cdot(1-\cos(\Delta\lambda))}{2}}\right)$ (18) $$

where $\Delta\varphi=\varphi_{2}-\varphi_{1}$, $\Delta\lambda=\lambda_{2}-\lambda_{1}$ and $r$ is taken to be $1$.

### E.1 Synthetic Regression

We utilise the GPyTorch software package (Gardner et al., ) for generating synthetic samples from a GP. As the number of datapoints in each sampled dataset is very large by GP standards ($1.1\times 10^{4}$), we approximate the SE kernel using structured kernel interpolation (SKI) (Wilson & Nickisch, ) with 100 grid points in each dimension. We use an observation noise of $\sigma_{n}=0.1$ for the smaller and larger lengthscale tasks. In Figure 14, we show an example dataset generated using a lengthscale of $0.1$ to demonstrate the complexity of these datasets. We were unable to compute ground truth log-likelihood values for these datasets without running into numerical issues.

Figure: (a) Ground-truth.
Refer to caption: extracted/5915947/figures/gp/ground-truth-gp-low-res.png

We train all models for $500,000$ iterations on $160,000$ pre-generated datasets using a batch size of eight. For all models, excluding the ConvCNP, we apply Fourier embeddings to each input dimension with $L=64$ wavelengths with $\lambda_{min}=0.01$ and $\lambda_{max}=12$. We found this to significantly improve the performance of all models. In Table 3, we provide test log-likelihood values for a number of gridded TNPs and baselines for both tasks. We observe that even when increasing the size of the baseline models they still underperform the smaller gridded TNPs. We include results for the Swin-TNP with full attention grid decoding (no NN-CA), which fail to model the more complex dataset when using the PT-GE.

**Table 3: Test log-likelihood ($\boldsymbol{\mathbf{\uparrow}}$) for the synthetic GP regression dataset. FPT: forward pass time for a batch size of eight in ms. Params: number of model parameters in units of M.**
| Model | Grid encoder | Grid size | $\ell=0.5$ ($\boldsymbol{\mathbf{\uparrow}}$) | $\ell=0.1$ ($\boldsymbol{\mathbf{\uparrow}}$) | FPT | Params |
| --- | --- | --- | --- | --- | --- | --- |
| CNP | - | - | $-0.406$ | $0.112$ | $9$ | $0.21$ |
| PT-TNP | - | $M=128$ | $0.819$ | $0.558$ | $53$ | $1.50$ |
| PT-TNP | - | $M=256$ | $0.819$ | $0.565$ | $74$ | $1.52$ |
| ConvCNP | SetConv | $32\times 32$ | $0.801$ | $0.536$ | $13$ | $2.11$ |
| ConvCNP | SetConv | $64\times 64$ | $0.830$ | $0.681$ | $93$ | $6.70$ |
| ViTNP | KI-GE | $32\times 32\rightarrow 16\times 16$ | $0.841$ | $0.722$ | $30$ | $1.16$ |
| ViTNP | PT-GE | $32\times 32\rightarrow 16\times 16$ | $0.841$ | $0.721$ | $32$ | $1.39$ |
| ViTNP | KI-GE | $16\times 16$ | $0.833$ | $0.711$ | $28$ | $1.09$ |
| ViTNP | PT-GE | $16\times 16$ | $0.840$ | $0.712$ | $29$ | $1.22$ |
| ViTNP | KI-GE | $64\times 64\rightarrow 32\times 32$ | $0.842$ | $0.728$ | $44$ | $1.16$ |
| ViTNP | PT-GE | $64\times 64\rightarrow 32\times 32$ | $0.836$ | $0.727$ | $56$ | $1.78$ |
| ViTNP | KI-GE | $32\times 32$ | $0.830$ | $0.725$ | $47$ | $1.09$ |
| ViTNP | PT-GE | $32\times 32$ | $0.837$ | $0.728$ | $53$ | $1.32$ |
| Swin-TNP | KI-GE | $32\times 32$ | $0.844$ | $0.723$ | $39$ | $1.09$ |
| Swin-TNP | PT-GE | $32\times 32$ | $0.844$ | $0.723$ | $42$ | $1.32$ |
| Swin-TNP | Avg-GE | $32\times 32$ | $0.840$ | $0.723$ | $34$ | $1.22$ |
| Swin-TNP | KI-GE | $64\times 64$ | $0.846$ | $0.728$ | $62$ | $1.09$ |
| Swin-TNP | PT-GE | $64\times 64$ | $\boldsymbol{\mathbf{0.847}}$ | $\boldsymbol{\mathbf{0.730}}$ | $69$ | $1.72$ |
| Swin-TNP | Avg-GE | $64\times 64$ | $0.845$ | $0.725$ | $58$ | $1.62$ |
| Swin-TNP (no NN-CA) | KI-GE | $32\times 32$ | $0.834$ | $0.716$ | $45$ | $1.09$ |
| Swin-TNP (no NN-CA) | PT-GE | $32\times 32$ | $0.837$ | $0.109$ | $48$ | $1.32$ |

##### ConvCNP

For the ConvCNP models, we use a regular CNN architecture with $C=128$ channels and five layers. We use a kernel size of five for the smaller ConvCNP ($32\times 32$ grid) and a kernel size of nine for the larger ConvCNP ($64\times 64$).

##### Swin-TNP

For the Swin-TNP models, we use a window size of $4\times 4$ for the smaller model ($32\times 32$ grid) and a window size of $8\times 8$ for the larger model ($64\times 64$ grid). The shift size is half the window size in each case. We also provide results when a simple average pooling is used for the grid encoder, which is similar to Xu et al. except that the pooling is performed in token space rather than on raw observations.

### E.2 Combining Weather Station Observations with Structured Reanalysis

Inspired by the real-life assimilation of 2m temperature (t2m), we use the ERA5 reanalysis dataset to extract skin temperature (skt) and 2m temperature (t2m) at a $0.25^{\circ}$ resolution (corresponding to a $721\times 1440$ grid). We then coarsen the skt grid to a $180\times 360$ grid, corresponding to $1^{\circ}$ in both the latitudinal and longitudinal directions. This implies that, because t2m lies on a finer grid, it essentially becomes an off-the-grid variable with respect to the coarsened grid on which skt lies. In order for the experimental setup to better reflect real-life assimilation conditions, we assume to only observe off-the-grid t2m values at real weather station locations(^12^1212More specifically, because the 2m temperature values come from the gridded ERA5 data, we only consider the nearest grid points to the true station locations as valid off-the-grid locations (i.e. if a station has coordinates at ($44.19^{\circ}$, $115.43^{\circ}$) latitude-longitude, we consider the grid point at ($44.25^{\circ}$, $115.5^{\circ}$))). In total, there are $9,957$ such weather station locations, extracted from the HadISD dataset (Dunn et al., ). We show their geographical location in Figure 15.

Figure: Figure 15: Station distribution within $1^{\circ}\times 1^{\circ}$ patches. The colour indicates the number of stations within each patch, clipped from a maximum value of 22 to 15. The distribution is far from uniform, with dense station areas in continents such as North America, Europe and parts of Asia, and a sparser distribution in Africa, South America, and in the oceans.
Refer to caption: extracted/5915947/figures/station_locations-low-res.png

For each task, we first randomly sample a time point, and then use the entire coarsened skt grid as the on-the-grid context data ($64,800$ points), as well as $N_{\text{off},c}$ off-the-grid t2m context points randomly sampled from the station locations. In the experiment from the main paper, $N_{\text{off},c}\sim\mathcal{U}_{[0,0.3]}$, but we also consider the case of richer off-the-grid context sets with $N_{\text{off},c}\sim\mathcal{U}_{[0.25,0.5]}$ in Table 6. The target locations are all the $9,957$ station locations.

We train all models for $300,000$ iterations on the hourly data between $2009-2017$ with a batch size of eight. Validation is performed on $2018$ and testing on $2019$. The test metrics are reported for $16,000$ data samples. Experiment specific architecture choices are described below.

##### Input embedding

We use spherical harmonic embeddings for the latitude / longitude values. These are not used in the ConvCNP model as the ConvCNP does not modify the inputs in order to maintain translation equivariance (in this case, with respect to the great-circle distance).

##### Grid sizes

For the main experiment (with results reported in Table 1), we chose a grid size of $64\times 128$ for the ConvCNP and Swin-TNP models, corresponding to a grid spacing of $2.8125^{\circ}$ in both the latitudinal and longitudinal directions.

In Table 6 we report results for a richer context set using a grid size of $128\times 256$ for the Swin-TNP models, corresponding to a grid spacing of $\approx 1.41^{\circ}$ in both the latitudinal and longitudinal directions. The results for the ConvCNP are for a grid size of $64\times 128$, to maintain a smaller gap in parameter count between models.

##### CNP

We use a different deepset for the on- and the off-the-grid data, and the mean as the permutation-invariant function to aggregate the context tokens, i.e. $\mathbf{z}_{c}=\frac{1}{N_{\text{off},c}}\sum_{n=1}^{N_{\text{off},c}}e_{\text
{off}}(\mathbf{x}_{\text{off},c,n},\mathbf{y}_{\text{off},c,n})+\frac{1}{N_{
\text{on},c}}\sum_{n=1}^{N_{\text{on},c}}e_{\text{on}}(\mathbf{x}_{\text{on},c
,n},\mathbf{y}_{\text{on},c,n})$

##### PT-TNP

We managed to use up to $M=256$ pseudo-tokens without running into memory issues. This shows that even if we only use two variables (one on- and one off-the-grid), PT-TNPs do not scale well to large data. We use a different encoder for the on- and off-the-grid data, before aggregating the two sets of tokens into a single context set.

##### ConvCNP

For the ConvCNP we use a grid of size $64\times 128$ for all experiments. We first separately encode both the on- and the off-the-grid to the specified grid size using the SetConv. We then concatenate the two and project them to a dimension of $C=128$ before passing through the U-Net (Ronneberger et al., ). The U-Net uses a kernel size of $k=9$ with a stride of one.

##### Swin-TNP

For the Swin-TNP models in the main experiment, we use a grid size of $64\times 128$, a window size of $4\times 4$ and a shift size of $2\times 2$.
For the experiment with richer off-the-grid context sets (i.e. between $0.25$ and $0.5$ of the off-the-grid data), we use a grid size of $128\times 256$, a window size of $8\times 8$ and a shift size of $4\times 4$.

#### E.2.1 Additional results for the main experiment

We provide in Figure 16 a comparison for an example dataset between the predictive errors (i.e. difference between predicted mean and ground truth) produced by three models: Swin-TNP with PT-GE, ConvCNP, and PT-TNP. The predictions are performed at all station locations. The stations included in the context set are indicated with a black dot. The figures show how Swin-TNP usually produces lower errors in comparison to the baselines, indicated through paler colours. Examples of regions where this is most prominent include central US, as well as southern Australia and southern Europe.

Figure: (a) Swin-TNP error.
Refer to caption: x5.png

##### Analysis of the predictive uncertainties

For the example dataset considered above, Figure 17 shows histograms of the normalised predictive errors, defined as the predictive errors divided by the predictive standard deviations. We compute the mean log-likelihoods under a standard normal distribution, and compare it to the reference negative entropy of the standard normal distribution of $-1.419$. This acts as an indicator of the accuracy of the predictive uncertainties outputted by the three models we consider: Swin-TNP (PT-GE, grid size $64\times 128$), ConvCNP (grid size $64\times 128$), and PT-TNP ($M=256$). For the dataset considered in Figure 17, we obtain $-1.439$ for Swin-TNP, $-1.490$ for ConvCNP, and $-1.380$ for PT-TNP, indicating that, out of the three models, Swin-TNP outputs the most accurate uncertainties.

Figure: (a) Swin-TNP histogram.
Refer to caption: extracted/5915947/figures/OOTG/histogram_swin.png

##### Analysis of grid size influence

We study to what extent increasing the grid size of the models, and hence their capacity, improves their predictive performance. We repeat the experiment for two models: Swin-TNP (with PT-GE) and ConvCNP with a grid size of $192\times 384$, corresponding to $0.9375^{\circ}$ in both latitudinal and longitudinal directions. For the Swin-TNP we use a window size of $8\times 8$, and a shift size of $4\times 4$. Due to time constraints, we only trained the Swin-TNP model for $280,000$ iterations instead of $300,000$. For the U-Net architecture within the ConvCNP we use a kernel size of nine. In the decoder of the bigger ConvCNP model, attention is performed over the nearest $9$ neighbours, whereas for the smaller ConvCNP we use full attention. This makes the FPT of the bigger model smaller that that of the $64\times 128$ model.

The results are shown in Table 4 (an represent a subset of the results shown in Table 1). In comparison to the Swin-TNP and ConvCNP models which use a grid size of $64\times 128$, both models improve significantly. However, the bigger ConvCNP still underperforms both the small and big variant of Swin-TNP, with a significant gap in both log-likelihood and RMSE.

**Table 4: Test log-likelihood ($\boldsymbol{\mathbf{\uparrow}}$) and RMSE ($\boldsymbol{\mathbf{\downarrow}}$) for the t2m station prediction experiment when varying grid size for two models. The standard errors of the log-likelihood are all below $0.004$, and of the RMSE below $0.005$. FPT: forward pass time for a batch size of eight in ms. Params: number of model parameters in units of M. Best results for each configuration (Swin-TNP / ConvCNP) are bolded.**
| Model | GE | Grid size | Log-lik. $\boldsymbol{\mathbf{\uparrow}}$ | RMSE $\boldsymbol{\mathbf{\downarrow}}$ | FPT | Params |
| --- | --- | --- | --- | --- | --- | --- |
| ConvCNP (no NN-CA) | SetConv | $64\times 128$ | $1.535$ | $1.252$ | 96 | $9.36$ |
| ConvCNP | SetConv | $192\times 384$ | $\boldsymbol{\mathbf{1.689}}$ | $\boldsymbol{\mathbf{1.166}}$ | 74 | $9.36$ |
| Swin-TNP | PT-GE | $64\times 128$ | $1.819$ | $1.006$ | 127 | $2.29$ |
| Swin-TNP | PT-GE | $192\times 384$ | $\boldsymbol{\mathbf{2.053}}$ | $\boldsymbol{\mathbf{0.873}}$ | $306$ | $10.67$ |

##### Analysis of influence of nearest-neighbour decoding

A final ablation we perform in this experiment is comparing the performance of two models with and without full attention at decoding time. More specifically, we compare Swin-TNP with PT-GE and KI-GE with a grid size of $64\times 128$, and either perform full attention in the grid decoder, or cross-attention over the $9$ nearest neighbours (NN-CA). For the variants with full attention, we evaluate the log-likelihood at $25\%$ randomly sampled station locations instead of all of them because of memory constraints. The results are shown in Table 5 (but are also presented in Table 1), and indicate that, not only does NN-CA offer a more scalable decoder attention mechanism, but it also leads to improved predictive performance when applied to spatio-temporal data. We hypothesise this is due to the inductive biases it introduces, which are appropriate for the strong spatio-temporal correlations present in the data we used.

**Table 5: Test log-likelihood ($\boldsymbol{\mathbf{\uparrow}}$) and RMSE ($\boldsymbol{\mathbf{\downarrow}}$) for the t2m station prediction experiment when varying the decoder attention mechanism—nearest-neighbour cross-attention and full attention (no NN-CA). The standard errors of the log-likelihood and RMSE are all below $0.003$. FPT: forward pass time for a batch size of eight in ms. Params: number of model parameters in units of M. Best results for each configuration (PT-GE / KI-GE) are bolded.**
| Model | GE | Grid size | Log-lik. $\boldsymbol{\mathbf{\uparrow}}$ | RMSE $\boldsymbol{\mathbf{\downarrow}}$ | FPT | Params |
| --- | --- | --- | --- | --- | --- | --- |
| Swin-TNP | KI-GE | $64\times 128$ | $\boldsymbol{\mathbf{1.683}}$ | $\boldsymbol{\mathbf{1.157}}$ | $121$ | $1.14$ |
| Swin-TNP (no NN-CA) | KI-GE | $64\times 128$ | $1.544$ | $1.436$ | $137$ | $1.14$ |
| Swin-TNP | PT-GE | $64\times 128$ | $\boldsymbol{\mathbf{1.819}}$ | $\boldsymbol{\mathbf{1.006}}$ | $127$ | $2.29$ |
| Swin-TNP (no NN-CA) | PT-GE | $64\times 128$ | $1.636$ | $1.273$ | $144$ | $2.29$ |

#### E.2.2 Additional Results for Richer Context Set

In order to investigate whether the model manages to learn meaningful relationships between the off- and on-the-grid data (t2m and skt) and to exploit the on-the-grid information(^13^1313This is achieved by comparing the performance of a model that is only given off-the-grid context information, with a similar model that is provided with both off- as well as on-the-grid data.), we also perform an experiment with richer context sets. More specifically, the number of off-the-grid context points is sampled according to $N_{\text{off},c}\sim\mathcal{U}_{[0.25,0.5]}$. The results are given in Table 6.

For the ConvCNP we evaluated two models—one with full decoder attention and one with nearest-neighbour cross-attention (NN-CA). Similarly to the previous section, we found that the latter has a better performance with a log-likelihood of $1.705$ (NN-CA) as opposed to $1.635$ (full attention). As such, Table 6 shows the results for the NN-CA ConvCNP model.

**Table 6: Test log-likelihood ($\boldsymbol{\mathbf{\uparrow}}$) and RMSE ($\boldsymbol{\mathbf{\downarrow}}$) for the t2m station prediction experiment with richer off-the-grid context information. The standard errors of both the log-likelihood and the RMSE are all below $0.003$. FPT: forward pass time for a batch size of eight in ms. Params: number of model parameters in units of M.**
| Model | GE | Grid size | Log-lik. ($\boldsymbol{\mathbf{\uparrow}}$) | RMSE ($\boldsymbol{\mathbf{\downarrow}}$) | FPT | Params |
| --- | --- | --- | --- | --- | --- | --- |
| CNP | - | - | $0.715$ | $3.056$ | $27$ | $0.34$ |
| PT-TNP | - | $M=256$ | $1.593$ | $1.403$ | $219$ | $1.57$ |
| ConvCNP | SetConv | $64\times 128$ | $1.705$ | $1.100$ | $21$ | $9.36$ |
| ViTNP | KI-GE | $48\times 96$ | $1.754$ | $1.112$ | $175$ | $1.14$ |
| ViTNP | PT-GE | $48\times 96$ | $1.988$ | $0.932$ | $188$ | $1.83$ |
| ViTNP | KI-GE | $144\times 288\rightarrow 48\times 96$ | $1.842$ | $1.046$ | $179$ | $1.29$ |
| ViTNP | PT-GE | $144\times 288\rightarrow 48\times 96$ | $2.242$ | $0.798$ | $221$ | $6.69$ |
| Swin-TNP | KI-GE | $128\times 256$ | $2.362$ | $0.758$ | $174$ | $1.14$ |
| Swin-TNP | PT-GE | $128\times 256$ | $\boldsymbol{\mathbf{2.446}}$ | $\boldsymbol{\mathbf{0.697}}$ | $208$ | $5.43$ |
| Swin-TNP (skt) | PT-GE | $128\times 256$ | $1.501$ | $1.266$ | $178$ | $5.43$ |
| Swin-TNP (t2m) | PT-GE | $128\times 256$ | $2.331$ | $0.909$ | $186$ | $5.38$ |

The results are consistent with the findings from the main experiment:

- •
The performances of the baselines (CNP, PT-TNP, and ConvCNP) are significantly worse than the gridded TNP variants considered.
- •
For each gridded TNP variant, the pseudo-token grid encoder (PT-GE) performs better than the kernel-interpolation one (KI-GE).
- •
The variants with a Swin-transformer backbone outperform the ones with a ViT-backbone, even when the latter has more parameters and a higher FPT.
- •
Among the ViT variants, the ones that employ patch encoding before projecting to a $48\times 96$ grid outperform the ones that directly encode to a $48\times 96$ grid.
- •
Performing nearest-neighbour cross-attention in the decoder as opposed to full attention leads to both computational speed-ups, as well as enhanced predictive performance.

In comparison to the main experiment, the gap in performance between Swin-TNP and Swin-TNP (t2m) is smaller—this is expected, given that the context already includes between $25\%$ and $50\%$ of the off-the-grid station locations. However, the gap is still significant, implying that Swin-TNP manages to leverage the on-the-grid data (skt) and exploits its relationship with the target t2m to improve its predictive performance.

### E.3 Combining Multiple Sources of Unstructured Wind Speed Observations

In this experiment, we consider modelling the eastward (u) and northward (v) components of wind speed at 700hPa, 850hPa and 1000hPa (surface level). These quantities are essential for understanding and simulating large-scale circulation in the atmosphere, for wind energy integration into power plants, or for private citizens and public administrations for safety planning in the case of hazardous situations (Lagomarsino-Oneto et al., ). We obtain each of the six modalities from the ERA5 reanalysis dataset (Hersbach et al., ), and construct datasets over a latitude / longitude range of $[25^{\circ},\ 49^{\circ}]$ / $[-125^{\circ},-66^{\circ}]$, which corresponds to the contiguous US, spanning four hours. We show plots of wind speeds at the three pressure levels for a single time step in Figure 18.

Figure: Figure 18: Wind-speed and direction at each of the three pressure levels, 700hPa, 850hPa and 1000hPa, over the contiguous US at 15:00 GMT, 08-06-1997. The colours correspond to the magnitude and direction of the wind speed.
Refer to caption: extracted/5915947/figures/wind-speed-us-matts-bday-low-res.png

Figure: (a) Swin-TNP error.
Refer to caption: extracted/5915947/figures/multi-modal/swintnp-24-60-1000hPa_pred_error_pc=0.05_ll=5.84_range=35.0-45.0_-85.0--70.0.png

Figure: (a) Swin-TNP normalised error.
Refer to caption: extracted/5915947/figures/multi-modal/swintnp-pt-24-60-1000hPa_pred_error_normalised_u_1000_ad=40.835_v_1000_ad=50.388_pc=0.05_ll=5.84_range=35.0-45.0_-85.0--70.0.png

For each task, we first sample a series of four consecutive time points. From this $4\times 96\times 236$ grid, we sample a proportion $p_{c}$ and $p_{t}$ of total points to form the context and target datasets, where $p_{c}\sim\mathcal{U}_{[0.05,0.25]}$ and $p_{t}=0.25$. All models are trained for $300,000$ iterations on hourly data between $2009-2017$ with a batch size of eight. Validation is performed on $2018$ and testing on $2019$. Experiment specific architecture choices are described below.

##### Input embedding

As the input contains both temporal and latitude / longitude information, we use both Fourier embeddings for time and spherical harmonic embeddings for the latitude / longitude values. These are not used in the ConvCNP as the ConvCNP does not modify the inputs to maintain translation equivariance (in this case, with respect to time and the great-circle distance).

##### Grid sizes

We chose a grid size of $4\times 24\times 60$ for the ConvCNP and Swin-TNP models, as this corresponds to a grid spacing of $1^{\circ}$ in the latitudinal direction and around $1^{\circ}$ in the longitudinal direction. For the ViTNP, we chose a grid size of $4\times 12\times 30$, corresponding to a grid spacing of $2^{\circ}$.

##### CNP

A different deepset is used for each modality, with the aggregated context token for each modality then summed together to form a single aggregated context token, i.e. $\mathbf{z}_{c}=\sum_{s=1}^{S}\frac{1}{N_{c,s}}\sum_{n=1}^{N_{c,s}}e_{s}(
\mathbf{x}_{c,n,s},\mathbf{y}_{c,n,s})$.

##### PT-TNP

For this experiment, we could only use $M=64$ pseudo-tokens for the PT-TNP without running into out-of-memory issues. This highlights a limitation in scaling PT-TNPs to large datasets. We note that there does exist work that remedies the poor memory scaling of PT-TNPs (Feng et al., ); however, this trades off against time complexity which itself is a bottleneck given the size of datasets we consider. A different encoder is used for each modality, before aggregating the tokens into a single context set of tokens.

##### ConvCNP

For the ConvCNP, we use a grid size of $4\times 24\times 60$. Each modality is first grid encoded separately using the SetConv, concatenated together and then to $C=128$ dimensions before passing through the U-Net.

##### Swin-TNP

For the Swin-TNP models, we use a grid size of $4\times 24\times 60$, a window size of $4\times 4\times 4$ and a shift size of $0\times 2\times 2$.

## Peer Reviews

confidence: 4
contribution: 2
initial_rating: 3
presentation: 3
questions: 1. How is the latent grid structured? That is, what resolution is chosen and how are the grid points distributed over the sphere, e.g., equirectangularly? Have you considered distortions towards polar regions in ERA5 data? I could find some information in Table 4 of the Appendix, but could not see for which tasks these grids were chosen (i.e., global vs local diagnosis).
2. What effect does the size of the attention circle (or diameter) in the decoder have on the performance? Depending on the target variable, this might be a crucial hyperparameter. However, increasing the attention radius increases the computation time. Can you provide some information on the relevance of the size of the attention circle for different target variables (e.g., t2m vs. wind) and what computational cost the size of the decoder's attention circle imposes?
soundness: 1
strengths: _Originality:_ The presented approach is original by combining established methods, such as Transformer Neural Processes, with efficient attention mechanisms from ViT and SwinTransformer, to tackle a highly relevant real-world problem. It would be very interesting to see, how the proposed method compares to ERA5 data and to understand the biases this approach imposes. Coming up with a data driven ERA5 alternative that mitigates numerical model biases would be of extraordinary value.

_Quality:_ The consideration of uncertainties in the output is a great design choice for the task of working with atmospheric states. Figures are mostly well designed and support the messages conveyed in the manuscript. Honestly, I do not see the value of Figure 1, in particular, since Figure 2 outlines the same pipeline both more detailed and concrete.

_Clarity:_ The manuscript is well organized, structured, and written. A rich appendix provides details about the model design, and a supplementary material with code is provided.
summary: A gridded transformer architecture is introduced that builds on neural processes. The model builds on three stages. First, unstructured data is encoded into a grid, second the grid is processed with an efficient ViT or SwinTransformer architecture, and third the grid is decoded to provide predictions at arbitrary positions beyond the grid. A comparison to other deep learning methods for processing unstructured grids demonstrates the benefits of the introduced method.
weaknesses: Overall, the paper touches on many interesting aspects, but misses a central message and problem that is tackled. Advertising the model more towards an alternative of traditional data assimilation techniques via numerical weather prediction models could be an interesting focus.

_Significance:_ It is hard to assess, whether the introduced method offers advantages over traditional data assimilation methods. To this end, the manuscript reads like introducing a neural network for data driven data aggregation (i.e., an alternative to ERA5), which potentially would provide substantial value. I would therefore suggest a focus shift from minimizing forecast error, such as done by all these deep learning weather prediction models, to generating a coherent state of the atmosphere, which is much more novel. Comparing to ERA5 would be highly appreciated, even though one cannot expect to outperform this dataset yet. Reporting first steps would be very insightful, though, and of large interest for the community.

_Quality:_
- **Qual1** Apart from the positive quality aspects mentioned under "Strengths," the manuscript can be improved in terms of supporting claims. That is, neural processes are introduced as flexible to downstream tasks, such as forecasting, data fusion, data interpolation and data assimilation (ll. 60--61). However, the model is only evaluated in terms of diagnosis error but not on the quality of other tasks.
- **Qual2** The use of different encoders for each input modality (l. 342) sounds inspired by ClimaX, which would be fair to be referenced here.
- **Qual3** Due to scaling reasons, the introduced method is not compared to its predecessors ANP and TNP (l. 366). However, it would be valuable to understand how the introduced method performs on a small dataset, where ANP and TNP are applicable, to understand whether the decign choices made in the new model result in performance decrease, and if so, to what extent.
- **Qual4** The manuscript advertises the transition from gridded reanalysis dataset, such as ERA5, from numerical weather prediction models. Yet, the model depends on ERA5 as additional gridded input beyond the heterogeneously distributed weather station data. It is unclear, how the gridded data and the weather station data is integrated and why the gridded data is necressary.
- **Qual5** GraphCast and GenCast follow a similar processing pipeline, where data points are aggregated into a grid, then processes, and subsequently decoded. A comparison to the aggregation approaches in these methods would be insightful. Otherwise it is hard to assess the quality of the introduced method over well established approaches.
- **Qual6** Figure 14 is very hard to interpret. A colourmap as well as a difference plot to the ground truth would help substantially to assess the quality of the different models. Similarly, it is difficult to extract the content of Figure 16, which is too dense. Zooming into certain regions could help.

_Clarity:_ Readers might benefit if the prediction task in Section 5.2 and Section 5.3 was introduced as diagnosis, i.e., diagnosing from some input stations in time $t$ to other input stations in the same time step $t$, where no forecasting is involved.

---

confidence: 3
contribution: 3
initial_rating: 5
presentation: 3
questions: See Weaknesses.
soundness: 3
strengths: The paper systematically introduces the G-TNPs model architecture, including the grid encoder, grid processor, and grid decoder, each accompanied by detailed mathematical descriptions and theoretical foundations.
summary: This paper introduces a Gridded Transformer Neural Process (G-TNP), which effectively maps unstructured spatiotemporal data onto a grid and incorporates an efficient attention mechanism, enhancing predictive performance on large-scale datasets.
weaknesses: 1. I'm curious, have you considered using large-scale meteorological models like Pangu as a baseline for comparison?
2. Although an efficient attention mechanism has been introduced, the model's complexity might still be high, especially when handling ultra-large-scale datasets. Could the authors provide comprehensive data to clarify this?

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: See above.
soundness: 3
strengths: Overall, I think the ability to blend off-grid and on-grid data of Gridded TNPs contributes well to the spatiotemporal modeling community. I like the simplicity of the proposed method and the clear writing of the paper. The experiments, even though are still at small scales, are reasonable for an academic setting, and sufficiently show the good performance of Gridded TNPs.
summary: The paper proposes Gridded Transformer Neural Processes (Gridded TNPs), an extension of TNPs that scales better to large-scale datasets. The idea is to use pseudo tokens and cross-attention to bring unstructured, point-based data to a grid before further processing by efficient on-grid architectures such as ViTs and Swin Transformers, and finally decoding to off-grid predictions via cross-attention. Different variants of Gridded TNPs were also introduced with different grid encoding methods and processor architectures. Gridded TNPs outperform the baselines on synthetic regression and station-based weather forecasting.
weaknesses: I do not see any critical weakness in the paper at the moment, but would like some clarifications and more analyses of the design choices in the paper. Specifically:
- What are the queries, keys, and values in the cross-attention operation in the grid encoder? What I understand from the paper is you use the learnable pseudo tokens as queries and embeddings of context points as keys and values. If this is true, do you use any positional encoding for the pseudo tokens?
- The paper said k=1 for the grid encoder is sufficient. Does this mean each pseudo token only attends to one context point that is closest to its grid location? What is the difference between this and simply picking the closest context point?
- In the grid decoder, where do the target tokens Z_t come from? Are they embeddings of the target inputs x's?
- I'd like to see more ablation studies about choosing k in both the encoder and decoder across different datasets and tasks.

----- POST REBUTTAL -----

The authors have addressed my main concerns so I'm raising my score to 8.

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: 1. GNNs are another type of commonly used approaches for handling unstructured data? Could the authors comment on them?
2. Compared to KI-GE, PT-GE appears to achieve better accuracy but requires a higher FPT. In the ViTNP cases in Table 1, the FPTs are 215 vs. 171 for PT-GE and KI-GE, respectively. Could the authors provide an explanation for this difference in time cost?
soundness: 3
strengths: The proposed gridded pseudo-token encoder and decoder with k-nearest-neighbors via cross-attention provided a flexible and elegant way to deal with both unstructured and structured data. The paper is well written, and the method is clearly presented.
summary: Transformers have shown promising results in weather forecasting but are limited to gridded structured data and neglect large amounts of unstructured observational data. 
Transformer Neural Processes (TNPs) are more flexible, capable of handling both structured and unstructured data, but they struggle with scalability. 
To address this issue, this paper introduces gridded pseudo-token TNPs.
Gridded TNPs consists of an encoder to convert the unstructured observation data to a structured grid via cross attention, a transformer (i.e., Vision Transformer [ViT] or Swin Transformer)-based processor to generate gridded token prediction, and a decoder with k-nearest-neighbor cross attention to map the predicted token back to the unstructured data at arbitrary locations.
The paper evaluates the methods on a few benchmarks, comparing the pseudo-token grid encoder to the kernel interpolation approaches and ViT to Swin Transformer.
weaknesses: 1. The experiment results are not analyzed in detail enough to substantiate the claims. The reviewer had a difficult time to link the very general discussions in section 5.4 to the experiments in sections 5.1-5.3. Adding a few sentences on the observations/conclusions of each experiment would be helpful.
2. The claim that introducing ViT and Swin Transformer into the transformer backbone as an efficient mechanism is a main contribution feels somewhat weak.
2. Typos: 
	1) In Fig. 2, the “{(x_{c,n},z_{c,n})}_n=1^N_c” should be “{(z_{c,n})}_n=1^N_c”?
	2) Line 344 on page 7, “e_s(x_{c,n,s},x_{c,n,s})” should be “e_s(x_{c,n,s},y_{c,n,s})”?

Meta-review: The paper presents transformer neural processes with specialized encoders and decoders to handle unstructured inputs and a attention-based processor to take advantage of the transformer scaling performance. They show experiments using ERA5 weather data by randomly sampling the data at sparse locations and using their model to predict temperature and wind at the same locations.

Strengths: A new method to deal with unstructured inputs is valuable to the community, well-organized, clear and detailed manuscript
Weaknesses: The biggest weakness is the limited evaluation of their model on tasks that require unstructured grids - the authors motivate the framework as general, broad, and flexible to a wide range of tasks in SciML or ML in general, but narrow the focus of their results on ERA5 and very simplified weather forecasting. Another smaller weakness is limited discussion of their model in comparison with many other ways to deal with unstructured data

Points for improvement:

Evaluation points:
The evaluation on subsampled ERA5 is positive and appreciated but the paper would have benefitted from stronger evaluations to motivate the contribution. The authors mention: Spatial interpolation e.g. downscaling,  Multi-resolution modelling, General regression tasks, Pattern completion, Feature prediction in arbitrary domains. Many of these have diverse datasets (or easily constructable), especially in PDE modeling (similar to weather forecasting). 

Discussion points: 
- discussions of GNNs (the request to compare to a GNN based model that includes GraphCast or graph neural operators is still valid since they should be able to process unstructured points, in principle - a discussion is necessary)
- discussions on ViT based models (transformers are designed to process arbitrary sequence of points - it is unclear from the paper why a ViT that processes a sequence of points cannot be used directly with a coordinate based position embedding or similar and needs to be atleast discussed). As noted in ClimaX referenced by one of the reviewers - "Moreover, since the standard ViT treats image modeling as pure sequence-to-sequence problems, it can perform tasks that some other variations cannot, such as learning from spatially incomplete data, where the input does not necessarily form a complete grid. This is useful in the regional forecasting task we consider in Section 4.2.2" - Pg 10 of https://arxiv.org/pdf/2301.10343.

The review points that requested the paper to be designed around data assimilation (replacement of ERA5, which is out-of-scope) or comparison with GenCast (retraining of the GenCast framework with using GraphCast to process an arbitrary mesh) were not used for the final decision.

Decision: Reject