# MAgNet: Mesh Agnostic Neural PDE Solver

Oussama Boussif
Mila - Québec AI Institute
DIRO, Université de Montréal
<EMAIL_ADDRESS>
&Dan Assouline
Mila - Québec AI Institute
DIRO, Université de Montréal
<EMAIL_ADDRESS>Loubna Benabbou
Université du Québec à Rimouski
<EMAIL_ADDRESS>&Yoshua Bengio
Mila - Québec AI Institute
DIRO, Université de Montréal
<EMAIL_ADDRESS>CIFAR Senior Fellow

###### Abstract

The computational complexity of classical numerical methods for solving
Partial Differential Equations (PDE) scales significantly as the resolution
increases. As an important example, climate predictions require fine spatio-
temporal resolutions to resolve all turbulent scales in the fluid simulations.
This makes the task of accurately resolving these scales computationally out
of reach even with modern supercomputers. As a result, current numerical
modelers solve PDEs on grids that are too coarse (3km to 200km on each side),
which hinders the accuracy and usefulness of the predictions. In this paper,
we leverage the recent advances in Implicit Neural Representations (INR) to
design a novel architecture that predicts the spatially continuous solution of
a PDE given a spatial position query. By augmenting coordinate-based
architectures with Graph Neural Networks (GNN), we enable zero-shot
generalization to new non-uniform meshes and long-term predictions up to 250
frames ahead that are physically consistent. Our Mesh Agnostic Neural PDE
Solver (MAgNet) is able to make accurate predictions across a variety of PDE
simulation datasets and compares favorably with existing baselines. Moreover,
MAgNet generalizes well to different meshes and resolutions up to four times
those trained on111Code and dataset can be found on:
https://github.com/jaggbow/magnet.

##  1 Introduction

Partial Differential Equations (PDEs) describe the continuous evolution of
multiple variables, e.g. over time and/or space. They arise everywhere in
physics, from quantum mechanics to heat transfer and have several engineering
applications in fluid and solid mechanics. However, most PDEs can’t be solved
analytically, so it is necessary to resort to numerical methods. Since the
introduction of computers, many numerical approximations were implemented, and
new fields emerged such as Computational Fluid Mechanics (CFD) (Richardson and
Lynch,, 2007). The most famous numerical approximation scheme is the Finite
Element Method (FEM) (Courant,, 1943; Hrennikoff,, 1941). In the FEM, the PDE
is discretized along with its domain, and the problem is transformed into
solving a set of matrix equations. However, the computational complexity
scales significantly with the resolution. For climate predictions, this number
can be quite significant if the desired error is to be reached, which renders
its use impractical.

In this paper, we propose to learn the continuous solutions for spatio-
temporal PDEs. Previous methods focused on either generating fixed resolution
predictions or generating arbitrary resolution solutions on a fixed grid (Li
et al.,, 2021; Wang et al.,, 2020). PDE models based on Multi-Layer
Perceptrons (MLPs) can generate solutions at any point of the domain
(Dissanayake and Phan-Thien,, 1994; Lagaris et al.,, 1998; Raissi et al.,
2017a, ). However, without imposing a physics-motivated loss that constrains
the predictions to follow the smoothness bias resulting from the PDE, MLPs
become less competitive than CNN-based approaches especially when the PDE
solutions have high-frequency information (Rahaman et al.,, 2018).

We leverage the recent advances in Implicit Neural Representations ((Tancik et
al.,, 2020), (Chen et al.,, 2020), (Jiang et al.,, 2020)) and propose a
general purpose model that can not only learn solutions to a PDE with a
resolution it was trained on, but it can also perform zero-shot super-
resolution on irregular meshes. The added advantage is that we propose a
general framework where we can make predictions given any spatial position
query for both grid-based architectures like CNNs and graph-based ones able to
handle sensors and predictions at arbitrary spatial positions.

##### Contributions

Our main contributions are in the context of machine learning for
approximately but efficiently solving PDEs and can be summarized as follows:

  * •

We propose a framework that enables grid-based and graph-based architectures
to generate continuous-space PDE solutions given a spatial query at any
position.

  * •

We show experimentally that this approach can generalize to resolutions up to
four times those seen during training in zero-shot super-resolution tasks.

##  2 Related Works

Current solvers can require a lot of computations to generate solutions on a
fine spatio-temporal grid. For example, climate predictions typically use
General Circulation Models (GCM) to make forecasts that span several decades
over the whole planet (Phillips,, 1956). These GCMs use PDEs to model the
climate in the atmosphere-ocean-land system and to solve these PDEs, classical
numerical solvers are used. However, the quality of predictions is
bottlenecked by the grid resolution that is in turn constrained by the
available amount of computing power. Deep learning has recently emerged as an
alternative to these classical solvers in hopes of generating data-driven
predictions faster and making approximations that do not just rely on lower
resolution grids but also on the statistical regularities that underlie the
family of PDEs being considered. Using deep learning also makes it possible to
combine the information in actual sensor data with the physical assumptions
embedded in the classical PDEs. All of this would enable practitioners to
increase the actual resolution further for the same computational budget,
which in turn improves the quality of the predictions.

Machine Learning for PDE solving:  Dissanayake and Phan-Thien, (1994)
published one of the first papers on PDE solving using neural networks. They
parameterized the solutions to the Poisson and heat transfer equations using
an MLP and studied the evolution of the error with the mesh size. Lagaris et
al., (1998) used MLPs for solving PDEs and ordinary differential equations.
They wrote the solution as a sum of two components where the first term
satisfies boundary conditions and is not learnable, and the second is
parameterized with an MLP and trained to satisfy the equations. In Raissi et
al., 2017a  the authors also parameterized the solution to a PDE using an MLP
that takes coordinates as input. With the help of automatic differentiation,
they calculate the PDE residual and use its MSE loss along with an MSE loss on
the boundary conditions. In follow-up work, Raissi et al., 2017b  also learn
the parameters of the PDE (e.g. Reynolds number for Navier-Stokes equations).

The recently introduced Neural Operators framework (Kovachki et al.,, 2021; Li
et al., 2020b, ; Li et al., 2020a, ) attempts to learn operators between
spaces of functions. Li et al., (2021) use "Fourier Layers" to learn the
solution to a PDE by framing the problem as learning an operator from the
space of initial conditions to the space of the PDE solutions. Their model can
learn the solution to PDEs that lie on a uniform grid while maintaining their
performance in the zero-shot super-resolution setting. In the same spirit,
Jiang et al., (2020) developed a model based on Implicit Neural
Representations called "MeshFreeFlowNet" where they upsample existing PDE
solutions to a higher resolution. They use 3D low-resolution space-time
tensors as inputs to a 3DUnet in order to generate a feature map. Next, some
points are sampled uniformly from the corresponding high-resolution tensors
and fed to an MLP called ImNet (Chen and Zhang,, 2018). They train their model
using a PDE residual loss and are able to predict the flow field at any
spatio-temporal coordinate. Their approach is closest to the one we propose
here. The main difference is that we perform super-resolution on the spatial
queries and forecast the solution to a PDE instead of only doing super-
resolution on the existing sequence.

Brandstetter et al., (2022) use the message-passing paradigm ((Gilmer et al.,,
2017), (Watters et al.,, 2017), (Sanchez-Gonzalez et al.,, 2020)) to solve 1D
PDEs. They are able to beat state-of-the-art Fourier Neural Operators (Li et
al.,, 2021) and classical WENO5 solvers while introducing the "pushforward
trick" that allows them to generate better long-term rollouts. Moreover, they
present an added advantage over existing methods since they can learn PDE
solutions at any mesh. However, they are not able to generalize to different
resolutions, which is a crucial capability of our method.

Most machine learning approaches require data from a simulator in order to
learn the required PDE solutions and that can be expensive depending on the
PDE and the resolution. Wandel et al., (2020) alleviate that requirement by
using a PDE loss.

Machine Learning for Turbulence Modeling:  Recent years have known a surge in
machine learning-based models for modeling turbulence. Since it is expensive
to resolve all relevant scales, some methods were developed that only solve
large scales explicitly and separately model sub-grid scales (SGS). Recently,
Novati et al., (2021) used multi-agent reinforcement learning to learn the
dissipation coefficient of the Smagorinsky SGS model (Smagorinsky,, 1963)
using as reward the recovery of the statistical properties of Direct Numerical
Simulations (DNS). Rasp et al., (2018) used MLPs to represent sub-grid
processes in clouds and replace previous parametrization models in a global
general circulation model. In the same fashion, Park and Choi, (2021) used
MLPs to learn DNS sub-grid scale (SGS) stresses using as input filtered flow
variables in a turbulent channel flow. Brenowitz and Bretherton, (2018) use
MLPs to predicts the apparent sources of heat and moisture using coarse-
grained data and use a multi-step loss to optimize their model.Wang et al.,
(2020) used one-layer CNNs to learn the spatial filter in LES methods and the
temporal filter in RANS as well as the turbulent terms. A UNet (Ronneberger et
al.,, 2015) is then used as a decoder to get the flow velocity. de Bezenac et
al., (2017) predict future frames by deforming the input sequence according to
the advection-diffusion equation and apply it to Sea-Surface Temperature
forecasting.

Stachenfeld et al., (2021) use the "encode-process-decode" (Sanchez-Gonzalez
et al.,, 2018, 2020) paradigm along with dilated convolutional networks to
capture turbulent dynamics seen in high-resolution solutions only by training
on low spatial and temporal resolutions. Their approach beats existing neural
PDE solvers in addition to the state-of-the-art Athena++ engine (Stone et
al.,, 2020). We take inspiration of this approach but replace the process
module by an interpolation module, to allow the model to capture spatial
correlations between known points and new query points.

##  3 Methodology

Figure 1: We illustrate the "Encode-Interpolate-Forecast" framework of MAgNet.
The parent mesh is fed to the encoder to generate the parent mesh embedding.
Next, we estimate the values at the spatial queries using the interpolation
module that uses features from both the parent mesh points and the parent mesh
embedding points closest to these queries. Finally, the parent mesh
observations and interpolated values at spatial queries are gathered as nodes
forming a new graph using nearest neighbors and the PDE solution is forecast
for all nodes (therefore all spatial locations) into the future using the
forecasting module.

We present the developed framework that leverages recent advances in Implicit
Neural Representations (INR) (Jiang et al.,, 2020; Sitzmann et al.,, 2020;
Chen et al.,, 2020; Tancik et al.,, 2020) and draws inspiration from mesh-free
methods for PDE solving. We first start by giving a mathematical definition of
a PDE. Next, we showcase the proposed "MAgNet" and derive two variants: A
grid-based architecture and a graph-based one.

###  3.1 Preliminaries

We define PDE as follows, using $D^{k}$ to denote $k$-th order derivatives:

######  Definition 3.1.

Evans, (2010) Let $U$ denote an open subset of $\mathbb{R}^{n}$ and $k\geq 1$
an integer. An expression of the form:

$\mathbf{\mathcal{L}}(D^{k}\mathbf{u}(x),D^{k-1}\mathbf{u}(x),\dots,\mathbf{u}(x),x)=\mathbf{0}\quad\forall
x\in U$ (1)

is called a $k$-th order system of PDEs, where
$\mathbf{\mathcal{L}}:\mathbb{R}^{mn^{k}}\times\mathbb{R}^{mn^{k-1}}\times\dots\times\mathbb{R}^{mn}\times\mathbb{R}^{m}\times
U\rightarrow\mathbb{R}^{m}$ is given and
$\mathbf{u}:U\rightarrow\mathbb{R}^{m},\mathbf{u}=(u^{1},\dots,u^{m})$ is the
unknown function to be characterized.

In this paper, we are interested in spatio-temporal PDEs. In this class of
PDEs, the domain is $U=[0,+\infty]\times\cal S$ (time $\times$ space) where
${\cal S}\subset\mathbb{R}^{n},\quad n\geq 1$ and, with $D^{k}$ indicating
differentiation w.r.t $x$, any such PDE can be formulated as:

$\begin{cases}\frac{\partial\mathbf{u}}{\partial
t}=\mathbf{\mathcal{L}}(D^{k}\mathbf{u}(x),\dots,\mathbf{u}(x),x,t)&\forall
t\geq 0,\forall x\in{\cal S}.\\\ u(0,x)=g(x)&\forall x\in\mathcal{S}\\\
\mathcal{B}u=0&\forall t\geq 0,\forall x\in\partial\mathcal{S}\end{cases}$ (2)

Where $\partial\mathcal{S}$ is the boundary of $\mathcal{S}$, $\mathcal{B}$ is
a non-linear operator enforcing boundary conditions on $u$ and
$g:\mathcal{S}\rightarrow\mathbb{R}^{m}$ represents the initial condition
constraints for the solution $u$.

Numerical PDE simulations have enjoyed a great body of innovations especially
where their use is paramount in industrial applications and research. Mesh-
based methods like the FEM numerically compute the PDE solution on a
predefined mesh. However, when there are regions in the PDE domain that
present large discontinuities, the mesh needs to be modified and provided with
many more points around that region in order to obtain acceptable
approximations. Mesh-based methods typically solve this problem by re-meshing
in what is called Adaptive Mesh Refinement (Berger and Oliger,, 1984; Berger
and Colella,, 1989). However, this process can be quite expensive, which is
why mesh-free methods have become an attractive option that goes around these
limitations.

###  3.2 MAgNet: Mesh-Agnostic Neural PDE Solver

####  3.2.1 "Encode-Interpolate-Forecast" framework

Let $\\{x_{1},x_{2},\dots,x_{T}\\}\in\mathbb{R}^{C\times N}$ denote a sequence
of $T$ frames that represents the ground-truth data coming from a PDE
simulator or real-world observations. $C$ denotes the number of physical
channels, that is the number of physical variables involved in the PDE and $N$
is the number of points in the mesh. These frames are defined on the same
mesh, that is the mesh does not change in time. We call that mesh the parent
mesh and denote its normalized coordinates of dimensionality $n$ by
$\\{p_{i}\\}_{1\leq i\leq N}\in[-1,1]^{n}$. Let $\\{c_{i}\\}_{1\leq i\leq
M}\in[-1,1]^{n}$ denote a set of $M$ coordinates representing the spatial
queries. The task is to predict the solution for subsequent time steps both
at: (i) all coordinates from the parent mesh $\\{p_{i}\\}_{1\leq i\leq N}$,
and (ii) coordinates from the spatial queries $\\{c_{i}\\}_{1\leq i\leq M}$.
At test time, the model can be queried at any spatially continuous coordinate
within the PDE domain to provide an estimate of the PDE solution at those
coordinates.

To perform the prediction, we first estimate the PDE solutions at the spatial
queries for the first $T$ frames and then use that to forecast the PDE
solutions at the subsequent timesteps at the query locations. We do this
through three stages (see Figure 1):

  1. 1.

Encoding: The encoder takes as input the given PDE solution
$\\{x_{t}\\}_{1\leq t\leq T}$ at each point of the parent mesh
$\\{p_{i}\\}_{1\leq i\leq N}$ and generates a state-representation of original
frames, which can be referred to as embeddings, and which we note
$\\{z_{t}\\}_{1\leq t\leq T}$. This representation will be used in the
interpolation step to find the PDE solution at the spatial queries
$\\{c_{i}\\}_{1\leq i\leq M}$. Note that in this encoding step, we can
generate one embedding for each frame such that we have $T$ embeddings or
summarize all the information in the $T$ frames into one embedding. We will
explain the methodology using $T$ embeddings, as it is easier to grasp the
time dimension in this formulation, but the implementation has been done using
a summarized single embedding, as mentioned in section 3.2.2. We also note
that the embedded mesh remains the same, i.e. we don’t change it by upsampling
or downsampling it.

  2. 2.

Interpolation: We follow the same approach as Jiang et al., (2020) and Chen et
al., (2020) by performing an interpolation in the feature space. Note that in
case we generate one representation that summarizes all $T$ frames into one,
then $z_{t}=z$ for $t=1,\dots,T$. Let $\\{t_{k}\\}_{1\leq k\leq T}$ denote the
timesteps at which the $x_{t}$ are generated.

For each spatial query $c_{i}$, let $\mathcal{N}(c_{i})$ denote the nearest
points in the parent mesh $p_{j}$. We generate an interpolation of the
features $z_{k}[c_{i}]$ at coordinates $c_{i}$ and at timestep $t_{k}$ as
follows:

$\forall k\in\\{1,T\\},\forall
i\in\\{1,M\\}:z_{k}[c_{i}]=\frac{\sum_{p_{j}\in\mathcal{N}(c_{i})}w_{j}g_{\theta}(x_{k}[p_{j}],z_{k}[p_{j}],c_{i}-p_{j},t_{k})}{\sum_{p_{j}\in\mathcal{N}(c_{i})}w_{j}}$
(3)

Where $z_{k}[p_{j}]$ and $x_{k}[p_{j}]$ denote the embedding and input frame
at position $p_{j}$ and time $t_{k}$ respectively. Moreover, $w_{j}$ are
interpolation weights and are positive and sum to one. Weights are chosen such
that points closer to the spatial query have a higher contribution to the
interpolated feature than points farther away from the spatial query. The
$g_{\theta}$ is an MLP. To get the PDE solution $x_{k}[c_{i}]$ at coordinate
$c_{i}$, we use a decoder $d_{\theta}$ which is an MLP here:
$x_{k}[c_{i}]=d_{\theta}(z_{k}[c_{i}])$. In practice, the number of neighbors
that we choose is $2^{n}$ where $n$ is the dimensionality of the coordinates.

  3. 3.

Forecasting: Now that we generated the PDE solution at the spatial queries
$c_{i}$ for all the past frames, we forecast the PDE solution at future time
points at both spatial queries and the parent mesh coordinates. Let
$\mathcal{G}$ denote the Nearest-Neighbors Graph (NNG) that has as nodes all
the $N$ locations in the parent mesh (at original coordinates
$\\{p_{i}\\}_{1\leq i\leq N}$ ) as well as all the $M$ query points (at
locations $\\{c_{i}\\}_{1\leq i\leq M}$), with edges that include only the
nearest neighbors of each node among the $N+M-1$ others. This corresponds to a
new mesh represented by the graph $\mathcal{G}$. Let
$\\{c^{\prime}_{i}\\}_{1\leq i\leq M+N}$ denote the corresponding new
coordinates. We generate the PDE solution for subsequent time steps on this
graph auto-regressively using a decoder $\Delta_{\theta}$ as follows:

$x_{k+1}[c^{\prime}_{i}]=x_{k}[c^{\prime}_{i}]+(t_{k+1}-t_{k})\Delta_{\theta}(x_{k}[c^{\prime}_{i}],\dots,x_{1}[c^{\prime}_{i}]),\quad
k=T,T+1,\dots$ (4)

We train MAgNet by using two losses:

  * •

Interpolation Loss: This loss makes sure that the interpolated points match
the ground-truth and is computed as follows:

$L_{\text{interpolation}}=\frac{\sum_{i=1}^{M}\sum_{k=1}^{T}||\hat{x}_{k}[c_{i}]-x_{k}[c_{i}]||_{1}}{T\times
M}$ (5)

Where $\hat{x}_{k}[c_{i}]$ denotes the interpolated values generated by the
model at the spatial queries.

  * •

Forecasting Loss: This loss makes sure that the model predictions into the
future are accurate. If $H$ is the horizon of the predictions, then we can
express the loss as follows:

$L_{\text{forecasting}}=\frac{\sum_{i=1}^{M+N}\sum_{k=1}^{H}||\hat{x}_{k+T}[c^{\prime}_{i}]-x_{k+T}[c^{\prime}_{i}]||_{1}}{H\times(M+N)}$
(6)

Where $\hat{x}_{k+T}[c^{\prime}_{i}]$ denotes the forecasted values generated
by the model at the graph $\mathcal{G}$ which combines both spatial queries
and the parent mesh.

The final loss is then expressed as:

$L=L_{\text{forecasting}}+L_{\text{interpolation}}.$

####  3.2.2 Implementation Details

In the previous section, we described the general MAgNet framework. In this
section, we present how we build the inputs to MAgNet as well as the
architectural choices for the encoding, interpolation and forecasting modules
and suggest two main architectures: MAgNet[CNN] and MAgNet[GNN].

Data pre-processing: We first consider a mesh that contains $N^{\prime}$
points ($N^{\prime}\geq N$). We randomly sample $N$ points from the mesh to
form the parent mesh. During training, $M$ spatial queries are randomly
sampled from the $N^{\prime}-N$ remaining points. We tried multiple values of
$M$ (that is the number of training spatial queries) to assess its impact on
the performance of the method within a sensitivity study presented in in
Section 4.3. The data pre-processing is illustrated in Figure 2

Figure 2: We illustrate the data pre-processing pipeline. We sample points
randomly from the starting mesh to form the parent mesh and the remaining
points form the spatial queries. Next, during training, we can sample from the
spatial queries and form what we call "training spatial queries". The
distinction is that the number of "training spatial queries" can be less than
the total number of spatial queries and we investigate the impact of this
number in Section 4.3.

MAgNet[CNN]: In this architecture, we follow Chen et al., (2020) and adopt the
EDSR architecture (Lim et al.,, 2017) as our CNN encoder. We concatenate all
frames $\\{x_{t}\\}_{1\leq t\leq T}$ in the channel dimension and feed that to
our encoder in order to generate a single representation $z$. For the
forecasting module, we use the same GNN as in (Sanchez-Gonzalez et al.,,
2020). A key advantage of this architecture is that it effectively turns
existing CNN architectures into mesh-agnostic ones by querying them at any
spatially continuous point of the PDE domain at test time.

MAgNet[GNN]: This model is similar to MAgNet[CNN] except that instead of using
a CNN as an encoder, we use a GNN: the same architecture as in the forecasting
module but each architecture having its separate set of parameters. This is
better suited for encoding frames with irregular meshes. Similarly to
MAgNet[CNN], we generate a single representation $z$ that summarizes all the
information from the frames $\\{x_{t}\\}_{1\leq t\leq T}$.

##  4 Results

In this section, we evaluate MAgNet’s performance against the following
baselines:

##### Fourier Neural Operators (FNO) (Li et al.,, 2021)

: Considered the state-of-the-art model in neural PDE solving, FNO casts the
problem of PDE solving as learning an operator from the space of initial
conditions to the space of the solutions. It is able to learn PDE solutions
that lie on a uniform grid and can do zero-shot super resolution.

##### Message-Passing Neural PDE Solvers (MPNN) (Brandstetter et al.,, 2022)

: Graph Neural Networks have been used to learn physical simulations with
great success (Sanchez-Gonzalez et al.,, 2020). Recently, they have been used
to learn solutions to PDEs (Brandstetter et al.,, 2022; Sanchez-Gonzalez et
al.,, 2020). MPNN-based GNNs coupled with an autoregressive strategy
demonstrate a superior performance to FNO and are able to make long rollouts
with the help of the "pushforward-trick" that only propagates gradient of the
last computed frame.

We evaluate all models on both 1D and 2D simulations (the datasets generated
from 1D and 2D PDEs are presented in Appendix A.4). All training sets in the
1D case contain 2048 simulations and test sets contain 128 simulations. For
the 2D case, training sets contain 1000 simulations and test sets contain 100
simulations. All models are evaluated using the Mean Absolute Error (MAE) on
the rolled out predictions averaged across time and space:

$MAE=\frac{\sum_{t=1}^{n_{t}-T}\sum_{i=1}^{N}|x_{t}[c_{i}]-\hat{x}_{t}[c_{i}]|}{(n_{t}-T)\times
N}.$ (7)

Where $n_{t}$ is the total number of frames and $T$ the number of frames input
to the models.

((a)) Zero-shot super-resolution performance on E1 test.

((b)) Zero-shot super-resolution performance on E2 test.

((c)) Zero-shot super-resolution performance on E3 test.

Figure 3: We present the models predictive performance on zero-shot super-
resolution, with a training spatial resolution of $n_{x}=50$ for the regular
grids. MAgNet[CNN] outperforms baselines on both E1 and E2 test sets but lags
behind FNO on E3. Error bars represent one standard deviation in both plots.

We train models for 250 epochs with early stopping with a patience of 40
epochs. See Appendix A.6 and A.7 for more implementation details.

We present a short summary of our results:

  * •

1D Case: Our model MAgNet is able to robustly generalize to unseen meshes in
the regular case as compared to the SOTA models FNO and MPNN. The performance
is even better in the irregular case. The details of the results are given in
the following section 4.1.

  * •

2D Case: Our model MAgNet is less competitive than MPNN when it comes to
generalizing to unseen meshes in the regular case. However, in the irregular
case, we are more competitive especially when trained on sparse meshes. The
results regarding the 2D PDEs are presented in section 4.2.

###  4.1 1D Case

For all the subsequent sections, $n_{x}$ and $n_{x}^{\prime}$ denote the
training and testing set’s resolutions respectively. The temporal resolution
$n_{t}=250$ remains unchanged for all experiments.

####  4.1.1 General Performance and zero-shot super resolution On Regular
Meshes

In this section, we compare MAgNet’s performance on all three datasets. All
models are trained on a resolution of $n_{x}=50$ and the PDE solutions lie on
a uniform grid. We test zero-shot super-resolution on
$n^{\prime}_{x}\in\\{40,50,100,200\\}$. Results are summarized in Figure 3 and
visualizations of the predictions can be found in Appendix A.2. MAgNet[CNN]
outperforms both baselines on both E1 and E2 datasets yet is slightly
outperformed by FNO on E3 (Figure 3(c)). Nonetheless, MAgNet[CNN]’s predictive
performance stays consistent up to $n^{\prime}_{x}=200$ while MPNN does not
generalize well to resolutions not seen during training.

####  4.1.2 General Performance and zero-shot super-resolution on Irregular
Meshes

In this section, we study how MAgNets compare against the other baselines when
it comes to making predictions on irregular meshes. In order to do so, we take
simulations from the uniform-mesh E1 dataset with a resolution of 100 do the
following:

Let $n_{x}\in\\{30,50,70\\}$. Then, for each simulation in the E1 dataset,
randomly sample the same subset of $n_{x}$ points: the mesh remains for each
single simulation in the E1 dataset.

The procedure is the same for the test set but we instead take the original E1
test set at a starting resolution of $200$ and generate four test sets with
irregular meshes for $n^{\prime}_{x}\in\\{40,50,100,200\\}$. This is different
from the test set of the previous section albeit considering the same
resolution since this one has irregular meshes. We summarize our findings in
Table 1.

MAgNet[GNN] performs better than MAgNet[CNN] on irregular meshes which is
expected since GNN encoders are better suited for this task. However,
surprisingly, even though we use a CNN encoder for MAgNet[CNN], the
performance seems to be better in most cases not only compared to FNO but also
MPNN which is a graph-based architecture. This effectively shows that MAgNet
can be used to turn existing CNN architectures into mesh-agnostic solvers.
This is particularly interesting for meteorological applications where one
needs to make predictions at the sub-grid level (at a specific coordinate)
while only having access to measurements on a grid.

| $\mathbf{n_{x}=30}$ | $\mathbf{n_{x}=50}$ | $\mathbf{n_{x}=70}$
---|---|---|---
Model | $\mathbf{n_{x}^{\prime}=40}$ | $\mathbf{n_{x}^{\prime}=50}$ | $\mathbf{n_{x}^{\prime}=100}$ | $\mathbf{n_{x}^{\prime}=200}$ | $\mathbf{n_{x}^{\prime}=40}$ | $\mathbf{n_{x}^{\prime}=50}$ | $\mathbf{n_{x}^{\prime}=100}$ | $\mathbf{n_{x}^{\prime}=200}$ | $\mathbf{n_{x}^{\prime}=40}$ | $\mathbf{n_{x}^{\prime}=50}$ | $\mathbf{n_{x}^{\prime}=100}$ | $\mathbf{n_{x}^{\prime}=200}$
FNO | 0.2784 | 0.2471 | 0.2574 | 0.2501 | 0.3797 | 0.3324 | 0.3841 | 0.3821 | 0.2798 | 0.2341 | 0.2533 | 0.2605
MAgNet[CNN] | 0.2081 | 0.1934 | 0.2063 | 0.2150 | 0.1869 | 0.1630 | 0.1599 | 0.1629 | 0.2237 | 0.1634 | 0.1385 | 0.1324
MPNN | 0.2602 | 0.1601 | 0.3451 | 0.3667 | 0.3027 | 0.2521 | 0.3226 | 0.3243 | 0.2685 | 0.1541 | 0.3403 | 0.3570
MAgNet[GNN] | 0.2422 | 0.2230 | 0.1938 | 0.1902 | 0.2302 | 0.1659 | 0.1590 | 0.1404 | 0.2400 | 0.1599 | 0.1398 | 0.1070

Table 1: We report the MAE per frame on the E1 dataset. We train all four
models on three different resolutions $n_{x}\in\\{30,50,70\\}$ and for each
training resolution, we evaluate zero-shot super-resolution on irregular
meshes for $n^{\prime}_{x}\in\\{40,50,100,200\\}$. We notice that even when we
use a CNN encoder, MAgNet not only performs better than the existing
baselines, but its performance stays consistent across different test
resolutions. MAgNet with a CNN encoder beats MPNN even when using an encoder
not suited for the task, which suggests MAgNet successfully turns existing CNN
architectures into mesh-agnostic ones.

###  4.2 2D case

In this section, we present results for the 2D PDE simulations. We use the
datasets B1 and B2 as our experimental testbeds (see appendix A.4). We use
$n_{train}$ to denote the train resolution and $n_{test}$ for the test
resolution. All models are fed a history of $T=10$ frames and are required to
generate a rollout of $n_{t}-T=40$ frames in the future.

####  4.2.1 General performance and zero-shot super resolution on regular
meshes

In this section, we compare MAgNet’s performance on both B1 and B2 datasets.
All models were trained on a resolution of $n_{train}=64^{2}$ and the PDE
solutions lie on a uniform grid. Zero-shot super resolution is tested on
$n_{test}\in\\{32^{2},64^{2},128^{2},256^{2}\\}$. We summarize our findings in
Figure 4. We notice that MAgNet[CNN] falls behind when it comes to making good
predictions and zero-shot super-resolution on both datasets while MAgNet[GNN]
and MPNN take the lead. Suprisingly, FNO suffers when it comes to generalizing
to unseen resolutions. Indeed, leveraging interactions between points after
the interpolation module allows MAgNet[GNN] not only to make good predictions
but also to generalize to denser and unseen resolutions.

Visualizations of the prediction for a sample for the B1 dataset for all test
resolutions can be found in Figures 10, 11, 12 and 13 respectively in the
appendix.

((a)) Results for the zero-shot super-resolution on the B1 dataset.

((b)) Results for the zero-shot super-resolution on the B2 dataset.

Figure 4: Results for the zero-shot super-resolution setting for models
trained in the regular setting on a resolution of $n_{train}=64^{2}$. While
FNO and MAgNet[CNN] fall behind, MPNN and MAgNet[GNN] take the lead by
leveraging the message-passing paradigm of Graph Neural Networks (Brandstetter
et al.,, 2022; Gilmer et al.,, 2017)

####  4.2.2 General performance and zero-shot super resolution on irregular
meshes

In this section, we present results for irregular meshes. We consider two
settings:

  * •

Uniform: $N$ nodes are randomly and uniformly sampled from a regular grid of
resolution $32^{2}$

  * •

Condensed: $N$ nodes are randomly sampled in a non-uniform way from a regular
grid of resolution $32^{2}$ following the distribution:
$p(x,y)\propto\exp\left(-8\left((x-0.25)^{2}+(y-0.25)^{2}\right)\right)$

We train MPNN and MAgNet[GNN] on four resolutions $N\in\\{64,128,256,512\\}$
where again $N$ is the number of nodes in the mesh for both Uniform and
Condensed settings. For each training setting, we evaluate both models on
regular grid PDE simulations from the B1 dataset. Figure 18 shows the mesh
nodes for the Uniform setting and Figure 19 is for the Condensed setting.

Our findings are summarized in Figure 5. We notice that our model MAgNet[GNN]
is especially more appealing in case we have fewer points to work with. In the
Condensed setting which is more realistic than the uniform one, MAgNet[GNN] is
even more competitive than MPNN. We note however, that while MAgNet[GNN]’s
performance is better for small test resolutions, it quickly deteriorates as
we test on higher resolutions which is a main limitation of our approach in
the 2D case.

((a)) Results for the zero-shot super-resolution on the irregular B1 dataset
for the Uniform setting.

((b)) Results for the zero-shot super-resolution on the irregular B1 dataset
for the Condensed setting.

Figure 5: Results for the zero-shot super-resolution setting for models
trained in the irregular setting on resolutions of $N\in\\{64,128,256,512\\}$.
MAgNet[GNN] is especially appealing in case we have fewer points and in case
the node distribution is non-uniform. However, the performance gap reduces as
both models see more points during training. We note however, that
MAgNet[GNN]’s performance deteriorates as we test on higher resolutions
compared to MPNN which is a main limitation of our approach in the 2D case.

###  4.3 Ablation study: Basic Interpolators vs Learned Interpolators

We investigate the contribution of the interpolation module to the general
predictive performance of MAgNet. We compare the MAgNet[CNN] architecture
against three ablated variants:

  * •

KNN: We use K-Nearest-Neighbors interpolation (Qi et al.,, 2017) on the
original frames directly to obtain the interpolated values at the spatial
queries.

  * •

Linear: We use Linear interpolation on the original frames directly to obtain
the interpolated values at the spatial queries.

  * •

Cubic: We use Cubic interpolation on the original frames directly to obtain
the interpolated values at the spatial queries.

Everything else is kept the same. The evaluation is done on the E1 dataset
with regular meshes and a resolution of $n_{x}=50$. Performance is tested on
E1 with regular meshes for test resolutions
$n^{\prime}_{x}\in\\{40,50,100,200\\}$. Results are summarized in Figure 6 and
additional ablation studies can be found in the appendix A.1.

Figure 6: We study the impact of having a learned interpolator (Ours) as
compared to existing interpolation schemes. Error bars represent one standard
deviation.

##  5 Limitations and Future Work

In this paper we introduced a novel framework that we call MAgNet for solving
PDEs on any mesh, possibly irregular. We proposed two variants of the
architecture which gave promising results on benchmark datasets. We were
effectively able to beat graph-based and grid-based architectures even when
using the CNN variant of the proposed framework, therefore suggesting a novel
way of adapting existing CNN architectures to make predictions on any mesh.
The main added value of the proposed method is its very good performance on
irregular meshes, including for super-resolution tasks, as it can be observed
for the presented 1D and 2D experiments, when compared to SOTA methods.
Notably, it seems to perform best when handed the smallest amount of data to
work with (with a 64x64 training resolution), even in the condensed non-
uniform case. This is a very desirable property for real-life data. A
limitation of our work however, is the significance of the learned
interpolator. Indeed, compared with a simple cubic interpolation, the approach
introduced here doesn’t seem to offer a significant advantage and we leave
improvement regarding this point for future work. Another improvement could be
seen in the forecasting module. For now, MAgNet forecasts using a first-order
explicit time-stepping scheme that is known to suffer from instability
problems in numerical PDE and ODE solvers. Learned solvers seem to somehow
circumvent this limitation even when using large time steps (Sanchez-Gonzalez
et al.,, 2020; Brandstetter et al.,, 2022; Stachenfeld et al.,, 2021). In a
future work, we wish to explore other time-stepping schemes such as the 4th
order Runge-Kutta method (Runge,, 1895; Kutta,, 1901) which is commonly used
for solving PDEs. Finally, we observed that our MAgNet model appears to
perform significantly better in a 1D than in a 2D setting (as shown in section
4.2). While it is retaining its superiority in an irregular mesh setting, it
appears a little less performant than the MPNN method, even with its GNN
variant. Further research on why the regular setting poses issues is thus also
reserved for future work.

## Acknowledgments and Disclosure of Funding

This work is financially supported by the government of Quebec and Samsung.
The authors would like to thank Shruti Mishra, Victor Schmidt, Dianbo Liu, and
Ayoub Ajarra for their fruitful discussions and useful insights.

## References

  * Bahdanau et al., (2015) Bahdanau, D., Cho, K., and Bengio, Y. (2015).  Neural machine translation by jointly learning to align and translate.  In Bengio, Y. and LeCun, Y., editors, 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings.
  * Bar-Sinai et al., (2018) Bar-Sinai, Y., Hoyer, S., Hickey, J., and Brenner, M. P. (2018).  Learning data driven discretizations for partial differential equations.
  * Berger and Colella, (1989) Berger, M. and Colella, P. (1989).  Local adaptive mesh refinement for shock hydrodynamics.  Journal of Computational Physics, 82(1):64–84.
  * Berger and Oliger, (1984) Berger, M. J. and Oliger, J. (1984).  Adaptive mesh refinement for hyperbolic partial differential equations.  Journal of Computational Physics, 53(3):484–512.
  * Brandstetter et al., (2022) Brandstetter, J., Worrall, D., and Welling, M. (2022).  Message passing neural pde solvers.
  * Brenowitz and Bretherton, (2018) Brenowitz, N. D. and Bretherton, C. S. (2018).  Prognostic validation of a neural network unified physics parameterization.  Geophysical Research Letters, 45(12):6289–6298.
  * Chen et al., (2020) Chen, Y., Liu, S., and Wang, X. (2020).  Learning continuous image representation with local implicit image function.
  * Chen and Zhang, (2018) Chen, Z. and Zhang, H. (2018).  Learning implicit fields for generative shape modeling.
  * Courant, (1943) Courant, R. (1943).  Variational methods for the solution of problems of equilibrium and vibrations.  Bulletin of the American Mathematical Society, 49(1):1–23.
  * de Bezenac et al., (2017) de Bezenac, E., Pajot, A., and Gallinari, P. (2017).  Deep learning for physical processes: Incorporating prior scientific knowledge.
  * Dissanayake and Phan-Thien, (1994) Dissanayake, M. W. M. G. and Phan-Thien, N. (1994).  Neural-network-based approximations for solving partial differential equations.  Communications in Numerical Methods in Engineering, 10(3):195–201.
  * Evans, (2010) Evans, L. (2010).  Partial Differential Equations.  American Mathematical Society.
  * Falcon et al., (2020) Falcon, W., Borovec, J., Wälchli, A., Eggert, N., Schock, J., Jordan, J., Skafte, N., Ir1dXD, Bereznyuk, V., Harris, E., Tullie Murrell, Yu, P., Præsius, S., Addair, T., Zhong, J., Lipin, D., Uchida, S., Shreyas Bapat, Schröter, H., Dayma, B., Karnachev, A., Akshay Kulkarni, Shunta Komatsu, Martin.B, Jean-Baptiste SCHIRATTI, Mary, H., Byrne, D., Cristobal Eyzaguirre, Cinjon, and Bakhtin, A. (2020).  Pytorchlightning/pytorch-lightning: 0.7.6 release.
  * Gilmer et al., (2017) Gilmer, J., Schoenholz, S. S., Riley, P. F., Vinyals, O., and Dahl, G. E. (2017).  Neural message passing for quantum chemistry.  In Proceedings of the 34th International Conference on Machine Learning - Volume 70, ICML’17, page 1263–1272. JMLR.org.
  * Hochreiter and Schmidhuber, (1997) Hochreiter, S. and Schmidhuber, J. (1997).  Long short-term memory.  Neural Computation, 9(8):1735–1780.
  * Hrennikoff, (1941) Hrennikoff, A. (1941).  Solution of problems of elasticity by the framework method.  Journal of Applied Mechanics, 8(4):A169–A175.
  * Jiang et al., (2020) Jiang, C. M., Esmaeilzadeh, S., Azizzadenesheli, K., Kashinath, K., Mustafa, M., Tchelepi, H. A., Marcus, P., Prabhat, and Anandkumar, A. (2020).  MeshfreeFlowNet: A Physics-Constrained Deep Continuous Space-Time Super-Resolution Framework.  IEEE Press.
  * Kingma and Ba, (2014) Kingma, D. P. and Ba, J. (2014).  Adam: A method for stochastic optimization.
  * Kovachki et al., (2021) Kovachki, N., Li, Z., Liu, B., Azizzadenesheli, K., Bhattacharya, K., Stuart, A., and Anandkumar, A. (2021).  Neural operator: Learning maps between function spaces.
  * Kutta, (1901) Kutta, W. (1901).  Beitrag zur näherungsweisen Integration totaler Differentialgleichungen.  Zeit. Math. Phys., 46:435–53.
  * Lagaris et al., (1998) Lagaris, I. E., Likas, A., and Fotiadis, D. I. (1998).  Artificial neural networks for solving ordinary and partial differential equations.  IEEE transactions on neural networks, 9 5:987–1000.
  * (22) Li, Z., Kovachki, N., Azizzadenesheli, K., Liu, B., Bhattacharya, K., Stuart, A., and Anandkumar, A. (2020a).  Neural operator: Graph kernel network for partial differential equations.
  * (23) Li, Z., Kovachki, N., Azizzadenesheli, K., Liu, B., Stuart, A., Bhattacharya, K., and Anandkumar, A. (2020b).  Multipole graph neural operator for parametric partial differential equations.  In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., and Lin, H., editors, Advances in Neural Information Processing Systems, volume 33, pages 6755–6766. Curran Associates, Inc.
  * Li et al., (2021) Li, Z.-Y., Kovachki, N. B., Azizzadenesheli, K., Liu, B., Bhattacharya, K., Stuart, A., and Anandkumar, A. (2021).  Fourier neural operator for parametric partial differential equations.  ArXiv, abs/2010.08895.
  * Lim et al., (2017) Lim, B., Son, S., Kim, H., Nah, S., and Lee, K. M. (2017).  Enhanced deep residual networks for single image super-resolution.
  * Novati et al., (2021) Novati, G., de Laroussilhe, H. L., and Koumoutsakos, P. (2021).  Automating turbulence modelling by multi-agent reinforcement learning.  Nature Machine Intelligence, 3(1):87–96.
  * Park and Choi, (2021) Park, J. and Choi, H. (2021).  Toward neural-network-based large eddy simulation: application to turbulent channel flow.  Journal of Fluid Mechanics, 914.
  * Paszke et al., (2019) Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Köpf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., and Chintala, S. (2019).  PyTorch: An Imperative Style, High-Performance Deep Learning Library.  Curran Associates Inc., Red Hook, NY, USA.
  * Phillips, (1956) Phillips, N. A. (1956).  The general circulation of the atmosphere: A numerical experiment.  Quarterly Journal of the Royal Meteorological Society, 82(352):123–164.
  * Qi et al., (2017) Qi, C. R., Yi, L., Su, H., and Guibas, L. J. (2017).  Pointnet++: Deep hierarchical feature learning on point sets in a metric space.
  * Rahaman et al., (2018) Rahaman, N., Baratin, A., Arpit, D., Draxler, F., Lin, M., Hamprecht, F. A., Bengio, Y., and Courville, A. (2018).  On the spectral bias of neural networks.
  * (32) Raissi, M., Perdikaris, P., and Karniadakis, G. E. (2017a).  Physics informed deep learning (part i): Data-driven solutions of nonlinear partial differential equations.
  * (33) Raissi, M., Perdikaris, P., and Karniadakis, G. E. (2017b).  Physics informed deep learning (part ii): Data-driven discovery of nonlinear partial differential equations.
  * Rasp et al., (2018) Rasp, S., Pritchard, M. S., and Gentine, P. (2018).  Deep learning to represent sub-grid processes in climate models.
  * Richardson and Lynch, (2007) Richardson, L. F. and Lynch, P. (2007).  Weather Prediction by Numerical Process.  Cambridge University Press.
  * Ronneberger et al., (2015) Ronneberger, O., Fischer, P., and Brox, T. (2015).  U-Net: Convolutional networks for biomedical image segmentation.  Medical Image Computing and Computer-Assisted Intervention – MICCAI 2015.
  * Runge, (1895) Runge, C. (1895).  Ueber die numerische auflösung von differentialgleichungen.  Mathematische Annalen, 46(2):167–178.
  * Sanchez-Gonzalez et al., (2020) Sanchez-Gonzalez, A., Godwin, J., Pfaff, T., Ying, R., Leskovec, J., and Battaglia, P. W. (2020).  Learning to simulate complex physics with graph networks.
  * Sanchez-Gonzalez et al., (2018) Sanchez-Gonzalez, A., Heess, N., Springenberg, J. T., Merel, J., Riedmiller, M., Hadsell, R., and Battaglia, P. (2018).  Graph networks as learnable physics engines for inference and control.
  * Sitzmann et al., (2020) Sitzmann, V., Martel, J. N. P., Bergman, A. W., Lindell, D. B., and Wetzstein, G. (2020).  Implicit neural representations with periodic activation functions.
  * Smagorinsky, (1963) Smagorinsky, J. (1963).  General circulation experiments with the primitive equations.  Monthly Weather Review, 91(3):99–164.
  * Stachenfeld et al., (2021) Stachenfeld, K., Fielding, D. B., Kochkov, D., Cranmer, M., Pfaff, T., Godwin, J., Cui, C., Ho, S., Battaglia, P., and Sanchez-Gonzalez, A. (2021).  Learned coarse models for efficient turbulence simulation.
  * Stone et al., (2020) Stone, J. M., Tomida, K., White, C. J., and Felker, K. G. (2020).  The athena++ adaptive mesh refinement framework: Design and magnetohydrodynamic solvers.  The Astrophysical Journal Supplement Series, 249(1):4.
  * Tancik et al., (2020) Tancik, M., Srinivasan, P. P., Mildenhall, B., Fridovich-Keil, S., Raghavan, N., Singhal, U., Ramamoorthi, R., Barron, J. T., and Ng, R. (2020).  Fourier features let networks learn high frequency functions in low dimensional domains.
  * Wandel et al., (2020) Wandel, N., Weinmann, M., and Klein, R. (2020).  Learning incompressible fluid dynamics from scratch – towards fast, differentiable fluid models that generalize.
  * Wang et al., (2020) Wang, R., Kashinath, K., Mustafa, M., Albert, A., and Yu, R. (2020).  Towards physics-informed deep learning for turbulent flow prediction.  In Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining. ACM.
  * Watters et al., (2017) Watters, N., Zoran, D., Weber, T., Battaglia, P., Pascanu, R., and Tacchetti, A. (2017).  Visual interaction networks: Learning a physics simulator from video.  In Guyon, I., Luxburg, U. V., Bengio, S., Wallach, H., Fergus, R., Vishwanathan, S., and Garnett, R., editors, Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc.

## Checklist

  1. 1.

For all authors…

    1. (a)

Do the main claims made in the abstract and introduction accurately reflect
the paper’s contributions and scope? [Yes]

    2. (b)

Did you describe the limitations of your work? [Yes]

    3. (c)

Did you discuss any potential negative societal impacts of your work? [No]

    4. (d)

Have you read the ethics review guidelines and ensured that your paper
conforms to them? [Yes]

  2. 2.

If you are including theoretical results…

    1. (a)

Did you state the full set of assumptions of all theoretical results? [N/A]

    2. (b)

Did you include complete proofs of all theoretical results? [N/A]

  3. 3.

If you ran experiments…

    1. (a)

Did you include the code, data, and instructions needed to reproduce the main
experimental results (either in the supplemental material or as a URL)? [Yes]
Link to the code and dataset is provided in the first footnot.

    2. (b)

Did you specify all the training details (e.g., data splits, hyperparameters,
how they were chosen)? [Yes] Training and implementation details can be found
in the Appendix

    3. (c)

Did you report error bars (e.g., with respect to the random seed after running
experiments multiple times)? [Yes] We ran models for 5 random seeds.

    4. (d)

Did you include the total amount of compute and the type of resources used
(e.g., type of GPUs, internal cluster, or cloud provider)? [Yes]

  4. 4.

If you are using existing assets (e.g., code, data, models) or
curating/releasing new assets…

    1. (a)

If your work uses existing assets, did you cite the creators? [Yes]

    2. (b)

Did you mention the license of the assets? [N/A]

    3. (c)

Did you include any new assets either in the supplemental material or as a
URL? [Yes] We provided a link to our repository that contains code and
datasets in the abstract.

    4. (d)

Did you discuss whether and how consent was obtained from people whose data
you’re using/curating? [N/A]

    5. (e)

Did you discuss whether the data you are using/curating contains personally
identifiable information or offensive content? [N/A]

  5. 5.

If you used crowdsourcing or conducted research with human subjects…

    1. (a)

Did you include the full text of instructions given to participants and
screenshots, if applicable? [N/A]

    2. (b)

Did you describe any potential participant risks, with links to Institutional
Review Board (IRB) approvals, if applicable? [N/A]

    3. (c)

Did you include the estimated hourly wage paid to participants and the total
amount spent on participant compensation? [N/A]

##  Appendix A Appendix

###  A.1 Additional Ablation and sensitivity studies

In this section, we study different architectural choices and the sensitivity
of key parameters in MAgNet.

##### Effect of modeling interactions between the parent mesh and the spatial
query

In section 3.2, we presented the developed method which can be used to
generate solutions at any spatial query through the "Encode-Interpolate-
Forecast" framework. This means that we are free to choose any architecture
for the three processes. In this section we investigate the choice of the
"Forecast" architecture on the predictive performance as well as zero-shot
super-resolution capabilities. We compare MAgNet[CNN] and a variant that uses
LSTM with attention ((Hochreiter and Schmidhuber,, 1997; Bahdanau et al.,,
2015) on the spatial queries only. Results are shown in Table 2. We see that
leveraging the interaction between the coordinates from the spatial query and
those in the parent mesh enables the model to give predictions consistent to
different resolutions as opposed to generating the solutions at these queries
with no interaction.

Table 2: Mean Absolute Error (MAE) reported for models trained on the E1 dataset with a resolution of $n_{x}=50$ on a uniform grid and tested on the same dataset for $n_{x}\in\\{40,50,100,200\\}$. We evaluate the effect of the Forecast module on the zero-shot super-resolution capabilities. The model without interaction contains an LSTM with attention (Bahdanau et al.,, 2015; Hochreiter and Schmidhuber,, 1997) for forecasting where spatial queries do not interact with the parent mesh. The model with interaction has a GNN that operates over the graph formed from the spatial-queries and the parent mesh (MAgNet[CNN]). Model | $n_{x}=40$ | $n_{x}=50$ | $n_{x}=100$ | $n_{x}=200$
---|---|---|---|---
Without Interaction | 0.1650 | 0.0815 | 0.2810 | 0.4139
With Interaction | 0.1079 | 0.1020 | 0.1142 | 0.1177

##### Impact of the number of point samples during training

We study the impact of the number of spatial queries used during training
($M$). We train MAgNet[GNN] on the E1 dataset with a resolution of $n_{x}=50$
on a uniform grid and test on the same resolution. Our findings are summarized
in Figure7. Increasing the number of spatial queries increases the predictive
performance as expected. Moreover, having many queries also decreases the
variance of the results. When we have fewer points, the random sampling can
cause some of these points to be in regions that decrease the loss faster than
other regions, hence, the model’s performance becomes sensitive to
randomization. However, this effect grows weaker as the number of queries
increases since they would uniformly cover more regions in the mesh.

Figure 7: We assess the impact of the number of spatial queries during
training. Error bars represent one standard deviation.

###  A.2 Visualizations: Models’ predictions in the 1D case

Figures 8 and 9 show the 1D models’ predictions on each of the test set
resolutions $n^{\prime}_{x}\in\\{40,50,100,200\\}$ on the E1 dataset.
MAgNet[CNN] predictions visually match the ground-truth’s while MPNN’s
prediction degrade as the predictions are advanced in time. The models shown
here are the ones trained on uniform meshes.

Figure 8: Vizualisation of the models’ 1D predictions on a simulation sample
from the E1 dataset. We present visualizations for each of the test
resolutions $n^{\prime}_{x}\in\\{40,50,100,200\\}$. The temporal resolution is
fixed at $n_{t}=250$. The x axis represents space and y axis represents time.
The arrow of time is from top to bottom. Figure 9: Vizualisation of the
models’ 1D predictions on a simulation sample from the E1 dataset. We present
visualizations for each of the test resolutions
$n^{\prime}_{x}\in\\{40,50,100,200\\}$. The temporal resolution is fixed at
$n_{t}=250$. The x axis represents space and y axis represents time. The arrow
of time is from top to bottom.

###  A.3 Visualizations: Models’ predictions in the 2D case

####  A.3.1 Regular case

Figures 10, 11, 12 and 13 show the 2D models’ predictions on each of the
methods on the B1 dataset, trained on regular meshes, for four different test
resolutions. The models are all trained on a 64$\times$64 grid. Visually MPNN
and MAgNet[GNN] predictions appear rather close, and both matching relatively
well the ground-truth’s.

Figure 10: Visualization of the models’ 2D predictions on a simulation sample
from the B1 dataset. The models are trained on a 64$\times$64 grid, and tested
on a 32$\times$32 resolution. The shown time steps are specified on top of the
figure. Figure 11: Visualization of the models’ 2D predictions on a simulation
sample from the B1 dataset. The models are trained on a 64$\times$64 grid, and
tested on a 64$\times$64 resolution. The shown time steps are specified on top
of the figure. Figure 12: Visualization of the models’ 2D predictions on a
simulation sample from the B1 dataset. The models are trained on a
64$\times$64 grid, and tested on a 128$\times$128 resolution. The shown time
steps are specified on top of the figure. Figure 13: Visualization of the
models’ 2D predictions on a simulation sample from the B1 dataset. The models
are trained on a 64$\times$64 grid, and tested on a 256$\times$256 resolution.
The shown time steps are specified on top of the figure.

####  A.3.2 Irregular case

Figures 14, 15, 16 and 17 show the 2D models’ predictions on the B1 dataset
for the Uniform setting on resolutions
$n_{test}\in\\{32^{2},64^{2},128^{2},256^{2}\\}$ respectively (see Section
4.2.2):

Figure 14: Visualization of the models’ 2D predictions on a simulation sample
from the B1 dataset, on an irregular mesh. The models are trained on a $N=512$
nodes, and tested on a 32$\times$32 resolution. The shown time steps are
specified on top of the figure. Figure 15: Visualization of the models’ 2D
predictions on a simulation sample from the B1 dataset, on an irregular mesh.
The models are trained on a $N=512$ nodes, and tested on a 64$\times$64
resolution. The shown time steps are specified on top of the figure. Figure
16: Visualization of the models’ 2D predictions on a simulation sample from
the B1 dataset, on an irregular mesh. The models are trained on a $N=512$
nodes, and tested on a 128$\times$128 resolution. The shown time steps are
specified on top of the figure. Figure 17: Visualization of the models’ 2D
predictions on a simulation sample from the B1 dataset, on an irregular mesh.
The models are trained on a $N=512$ nodes, and tested on a 256$\times$256
resolution. The shown time steps are specified on top of the figure.

###  A.4 PDE Datasets

##### 1D PDE Simulations

For the 1D case, We use three of MPNN’s PDE simulations (Brandstetter et al.,,
2022) as our experimental testbed. In the same fashion, we are interested in
the following family of PDEs:

$\begin{cases}[\partial_{t}u+\partial_{x}(\alpha
u^{2}-\beta\partial_{x}u+\gamma\partial_{xx}u)](t,x)=\delta(t,x)&\\\
u(0,x)=\delta(0,x)&\\\ \delta(t,x)=\sum_{j=1}^{J}A_{j}\sin(\omega_{j}t+2\pi
l_{j}x/L+\phi_{j})\end{cases}$ (8)

Where $J=5$, $L=16$ and coefficients sampled uniformly in
$A_{j}\in[-0.5,0.5]$, $\omega_{j}\in[-.4,-0.4]$, $l_{j}\in\\{1,2,3\\}$,
$\phi_{j}\in[0,2\pi)$ and periodic boundary conditions following Brandstetter
et al., (2022); Bar-Sinai et al., (2018). For the 1D case, the temporal
resolution is set to $n_{t}=250$ for the entire study. During testing, all the
models are fed a history of $T=25$ frames and produce a rollout of
$n_{t}-T=225$ frames in the future. Here, we present the two datasets we work
with:

  * •

E1: Burgers equation without diffusion $\alpha=1,\beta=0,\gamma=0$

  * •

E2: Burgers equation with variable diffusion $\alpha=1,\beta=\eta,\gamma=0$
where $\eta\in[0,0.2]$

  * •

E3: Mixed scenario where $\alpha\in[0,3],\beta\in[0,0.4]$ and
$\gamma\in[0,1]$.

##### 2D PDE Simulations

For the 2D case, we use the Burgers equation as our experimental testbed. This
dataset was generated using Phiflow 222https://github.com/tum-pbs/PhiFlow:

$\begin{cases}[\partial_{t}u+u\nabla u](t,x,y)=\beta[\Delta u](t,x,y)&\\\
u(0,x,y)=f(x,y)&\\\ f(x,y)=\sum_{j=1}^{J}A_{j}\sin(2\pi
l_{j}^{x}x/L+\phi_{j}^{x})\cos(2\pi l_{j}^{y}y/L+\phi_{j}^{y})\end{cases}$ (9)

Where $J=5$, $L=64$ and coefficients sampled uniformly in
$A_{j}\in[-0.5,0.5]$, $l_{j}^{x},l_{j}^{y}\in\\{1,2,3\\}$,
$\phi_{j}^{x},\phi_{j}^{y}\in[0,2\pi)$ and periodic boundary conditions. For
the 2D case, the temporal resolution is set to $n_{t}=50$ for the entire
study. During testing, all the models are fed a history of $T=10$ frames and
produce a rollout of $n_{t}-T=40$ frames in the future. Here, we present the
three datasets we work with:

  * •

B1: Burgers equation without diffusion $\beta=0$

  * •

B2: Burgers equation with variable diffusion $\beta\in(0,0.2]$

###  A.5 Visualization of the irregular meshes in the 2D case

Figure 18: Visualization of an irregular mesh, uniformly distributed. Figure
19: Visualization of an irregular mesh, non-uniformly distributed but
following a bi-variate normal distribution.

###  A.6 Training details

We train all models for 250 epochs and early stopping with a patience of 40
epochs. All models are trained using Adam Optimizer (Kingma and Ba,, 2014) and
the StepLR learning scheduler (Paszke et al.,, 2019) which decays the learning
rate by a factor $k$ every $N_{steps}$ epochs. All models were trained on 5
random seeds $\in\\{5,10,21,42,2022\\}$. We use Pytorch (Paszke et al.,, 2019)
and Pytorch-Lightning (Falcon et al.,, 2020) for our implementations and we
summarize all model’s training hyper-parameters in Table 3.

Table 3: Training hyperparameters for FNO, MPNN, MAgNet[CNN] and MAgNet[GNN] Parameters | FNO | MPNN | MAgNet[CNN] | MAgNet[GNN]
---|---|---|---|---
Learning Rate | 0.001 | 0.001 | 0.001 | 0.001
Weight Decay | 0 | 0 | 0 | 0
$k$ | 0.3 | 0.3 | 0.3 | 0.3
$N_{steps}$ | 50 | 50 | 40 | 50
GPU | RTX8000 | RTX8000 | RTX8000 | RTX8000
Number of GPUs | 1 | 1 | 2 | 2
Training duration (hours) | 1 | 1 | 5 | 2.5

###  A.7 Architectural details

####  A.7.1 MAgNet[CNN]

##### Encoder Architecture

We adapt the original EDSR (Lim et al.,, 2017) architecture to work on 1D
signals instead of 2D and use 4 residual blocks with a hidden dimension of
$128$.

##### Interpolation Module

We use a $4$ layers MLP with a hidden size of $64$ followed by Layernorm for
$g_{\theta}$. For $d_{\theta}$, we use a $4$ layer MLP with a hidden size of
$64$

##### Forecasting Module

We use the same architecture as in Sanchez-Gonzalez et al., (2020). The
encoder module uses a $4$ layer MLP with a hidden size of $64$ and the latent
dimension to $32$. We use $5$ message-passing steps (with the same parameters
for the MLP), the decoder also has a $4$ layer MLP with hidden size of $64$.

####  A.7.2 MAgNet[GNN]

##### Encoder Architecture

We use the same architecture as in Sanchez-Gonzalez et al., (2020) but only
keep the encoder and processor. We also use 5 message-passing steps for the
processor and use $4$ layer MLP with hidden size of $128$ and a latent
dimension of $128$.

##### Interpolation Module

We use a $4$ layers MLP with a hidden size of $128$ followed by Layernorm for
$g_{\theta}$. For $d_{\theta}$, we use a $4$ layer MLP with a hidden size of
$128$

##### Forecasting Module

We use the same architecture as in Sanchez-Gonzalez et al., (2020). The
encoder module uses a $4$ layer MLP with a hidden size of $64$ and the latent
dimension set to $128$. We use $5$ message-passing steps for the processor
(with the same parameters for the MLP), the decoder also has a $4$ layer MLP
with hidden size of $128$.

####  A.7.3 MPNN

We use the same hyperparameters as in Brandstetter et al., (2022).

####  A.7.4 FNO

We use $5$ Fourier Layers and $12$ modes in Fourier space. We use a hidden
channel size of $256$.


††thanks:  This work was supported by the U.S. Department of Energy, Office of
Basic Energy Sciences under contract DE-AC02-06CH11357.

# Precise Orbital Tracking of an Asteroid with a Phased Array of Radio
Transponders

Bernhard W. Adams  Argonne National Laboratory, 9700 S. Cass Ave., Argonne, IL
60439, USA<EMAIL_ADDRESS>
###### Abstract

Deflecting an asteroid from an Earth impact trajectory requires only small
velocity changes, typically of the order of microns per second, if done many
years ahead of time. For this, a highly precise method of determining the
need, magnitude, and direction of a deflection is required. Although the
required precision can be achieved by much less accurate extended
observations, an intrinsic resolution of $\mu$m/s permits the live monitoring
of nongravitational orbital perturbations (Yarkovsky effect), and of a
deflection effort itself. Here, it is proposed to deploy on the asteroid’s
surface multiple radio units to form a phased array capable of measuring
radial velocities relative to Earth to about 1 $\mu$m/s and ranges to 5 m. The
same technology can also be used for scientific applications such as very-long
baseline radio astronomy, milli-Hertz gravitational wave detection, or mapping
of the solar wind.

##  I Introduction

Asteroid impacts on Earth are infrequent, but potentially devastating events.
One example is (99942) Apophis, a 300-meter rock, which will miss Earth in
2029 by less than the orbital height of geostationary satellites. If, in this
encounter, its center of gravity passes through the so-called keyhole – an
imaginary, 600-meter wide region in space near Earth – then gravitational
deflection will set it up for an impact in 2036 with an energy release similar
to the explosion of Krakatoa in 1883. Making the asteroid miss the keyhole
requires much less effort than avoiding the entire planet once the keyhole has
been traversed: Up to three years before 2029, a velocity change $\Delta
v<1\mu$m/s tangential to the orbit is sufficient to avoid the keyhole san ;
Schweickart et al. (2006), but a $\Delta v$ of the order of cm/s is necessary
for steering away from the entire planet after a keyhole passage. Therefore, a
determination whether a deflection will be necessary should be done well
before 2029. To obviate the need for an unnecessarily large deflection, the
orbit should be measured to a precision commensurate with the required
velocity change, i.e., 1 $\mu$m/s. Passive radar can resolve velocities to a
few cm/s and ranges to tens of meters, and only at distances much less than
the orbital diameter of Earth and Apophis. A much higher accuracy can be
achieved by use of a radio transponder traveling with or on the asteroid.

Placement and operation of a radio transponder on an asteroid is challenging
in several ways, including (i) deliverable payload, (ii) power supply, (iii)
unknown surface topography and composition, (iv) landing without bouncing off
faster than the escape velocity of a few cm/s, (v) maintaining a directional
link to Earth while the asteroid is rotating, (vi) degradation of radio
signals due to interplanetary-plasma scintillations, and (vii) the space
environment. The problem of landing would be irrelevant, and the power supply
would be simpler with a transponder aboard a spacecraft orbiting/accompanying
the asteroid, but transponders fixed to the asteroid surface can reach a much
higher accuracy (see below). Therefore, the present proposal will concentrate
on the latter approach.

##  II Design Goals

First, to the accuracy: Doppler tracking of planetary probes has shown
performance at the stipulated level, most recently 0.5 $\mu$m/s for the
gravitational-wave experiment Asmar et al. (2005) on Cassini currently
orbiting Saturn. Orbital refinement techniques achieve this with much less
precise individual measurements over an orbital period, but a high-precision
technology permits live monitoring of a deflection maneuver, as well as of the
Yarkovsky Chesley et al. (2003) effect (orbital perturbations due to radiation
pressure). This would also be useful for using the Yarkovsky itself for
deflection Spitale (2002). In addition to Doppler tracking, this proposal also
aims for ranging at 5 m resolution, corresponding to 1 $\mu$m/s over two
months, about 20% of Apophis’ orbital period. These Doppler and ranging
accuracies are stipulated for most of the asteroid’s orbit, except where the
line of sight passes within 100 solar radii from the sun, and radio-signal
degradation due to scintillations (see below) becomes too severe. This makes
the maximum tracking distance 1.82 astronomical units (A.U.), i.e., 2.72
$\cdot 10^{8}$ km.

A radio transponder could be placed on the surface or aboard a companion
spacecraft orbiting or tracking the asteroid. The latter option can be
realized with proven spacecraft technology, but it cannot reach $\mu$m/s
accuracy for the asteroid’s center of mass in a single measurement: due to the
irregular shape and rotation of the asteroid, orbits are generally nonperiodic
and unstable. Therefore, long averaging is necessary without orbit correction
manouvers, but these would be required to prevent crashing or ejection. Active
locking of a companion spacecraft to surface features of the asteroid using
laser or radar ranging is also unlikely to be accurate to the $\mu$m/s level.
Therefore, instantaneous high-precision measurements for live monitoring of
deflections, or of the Yarkovsky effect, require a transponder on the
asteroid’s surface.

##  III Details of the Design

For reasonable power consumption, a directional antenna is required that
compensates for the asteroid’s rotation. Using a tracking dish antenna on the
low-gravity surface seems challenging. A much better approach is found in
phased-array technology, where the phases of received and transmitted signals
in a multitude of dipole antennae are controlled electronically without any
mechanical motion for directionality through constructive interference. The
phased array discussed here Adams (2007) consists of 475 radio units grouped
in 19 nodes of 25 each (Fig. 1). These numbers are somewhat arbitrary, and are
meant only as an example.

Figure 1: Left: The large web connecting 19 nodes by optical fibers and thin
electrical wires. Right: One node containing 25 transponders, connected to
each other by optical fibers and thin metal wires attached to a thin plastic
foil. As shown, the total length of wire connecting the nodes is 11.4 km.
Lower right: cross section of a wire with carbon fibers, optical fibers, and
metal filaments.

The nodes and the radio units within each node are linked together by a net of
thin threads containing carbon fibers for mechanical strength, glass fibers
for the distribution of a precise timing signal (for phasing) and digital
communication, and metal wires for electric-power sharing. Each radio unit has
24 antennae for a total of 11400 for each of the frequencies used (see below).
Each node contains a highly stable oscillator, one of which generates a
reference frequency that is distributed over the glass-fiber network. Shortly
before arrival, the transfer spacecraft aims for a collision with the asteroid
at a velocity of about 0.1 m/s (roughly the escape velocity from the
asteroid’s surface). As shown in Fig. 2, it then sequentially ejects 18
containers

Figure 2: Deployment of the coarse web. A space vehicle (1) carries the array
to the asteroid. It is comprised of three sections, each containing six
slices. These are ejected in sequence (2-3-4) to deploy the web. After
ejection of all slices, the entire web is made to rotate slowly to keep it
stretched out (see text). Finally, each of the slices deploys a sub-web, as
shown in fig. 3.

for the nodes, which are connected by carbon/metal/glass-fiber threads. To
keep the fully deployed net spread out over about 800 meters, it is then spun
up by gas jets to about one revolution per hour. Next, each container
(including one in the remainder of the spacecraft) releases a thin plastic
foil with attached radio units and inflates it to 50 m diameter by injecting
gas into thin tubes (Fig. 3).

Figure 3: Ejection and unfolding of the sub-webs: 1: container opens, 2:
folded plastic foil is ejected, 3: foil unfolds as gas is blown into tubes
formed between two layers of the foil.

Continuing on its collision course, the assembly wraps itself around the
asteroid. After several slow-motion bounces that are restrained by the net,
everything settles down in fixed random locations. Some damage is tolerable as
it only leads to a slightly degraded performance (see below).

Next, Earth sends a radio signal sufficiently strong for reception without
phasing (see below). As the asteroid rotates, the array compares the phases of
this signal in each antenna to the local reference, and thus acquires the
information necessary for phased-array operation at each of the frequencies
used. This information is continually refined as the asteroid and Earth
progress along their orbits.

The Doppler-tracking resolution is degraded by scintillations due to
fluctuations of the radio refractive index n of the turbulent interplanetary
plasma (solar wind). Extrapolating spacecraft-communication data Woo (1978)
for a frequency of $\nu=2.3$ GHz from 20 to 100 solar radii, and applying
scaling Ho et al. (2002) with $\nu^{-1.2}$, the spectral broadening due to
scintillations of an 8-GHz signal at 100 solar radii is about 0.01 Hz. Ways to
mitigate them are discussed in the literature on spacecraft-based
gravitational-wave experiments Armstrong et al. (2003); Armstrong (2006),
using the dispersion $n-1\sim n^{-2}$ of $n$ with multiple frequencies $\nu$
to measure, and then correct for the fluctuations. Here, frequencies of about
8 GHz and 16 GHz will be used (slightly offset in up- and downlink).

Doppler shifts are detected through phase shifts between ground and remote
clocks, and a radial-velocity resolution of $\Delta v_{r}=1\mu$m/s requires
the reference clocks to be stable to $\Delta v_{r}/c=3\cdot 10^{-15}$, where
$c$ is the speed of light. With a radio frequency of $\nu=8\mbox{ GHz}$ and a
phase sensitivity of $s=2\pi/100$ (3.5 degrees), this stability must be
maintained over the time $\tau=(c/\Delta v_{r})s/(2\pi\nu)=375$ s. This
exceeds the state of the art of 1 part in $10^{-13}$ over 1000 s in spaceborne
ultra-stable oscillators (USOs), but can be achieved with a coherent link
Howard et al. (1974); Woo et al. (1976), where the uplink signal is used to
phase-lock the space-borne oscillator from which the return signal is derived.
Recently, the modified technique of noncoherent Doppler tracking DeBoy et al.
(2005) was used, where the spacecraft has a freely-running oscillator and
transmits its phase slippage relative to the signal received from Earth
through a digital communication channel. This latter technique is proposed
here Adams (2007) to make it easier to accommodate radio signals at two
frequencies that fluctuate differently as they traverse the solar wind. With
the above oscillator stability of $10^{-13}$, the 1-$\sigma$ walk-off of an
8-GHz signal is $2\pi$ in 1250 s. Thus, to resolve the phase to $s=2\pi/100$,
about 8 bits of digital data need to be transmitted every 1250 s for
1-$\sigma$ certainty, and 10 bits for 4-$\sigma$.

For ranging, Earth and the array transmit radio pulses to each other with a
pseudo-random number (PRN) phase modulation similarly to GPS, and cross-
correlate the received signals with the local ones. The array transmits the
correlation result through the digital communication channel. With Gaussian
Minimum-Shift Keying (GMSK) and a bit-time-bandwidth product of BT = 0.3
(standard in cell phonesAgilent Technologies (2005)), the signal bandwidth for
a 5-meter resolution is about 20 MHz Adams (2007). Here, $\nu=2.25$ GHz is
chosen for the uplink, and $\nu=2.1$ GHz for downlink ranging. The pulse
duration must cover the initial ranging uncertainty of about 1 km, that is 3
$\mu$s, but shorter PRN sequences may be used once the range is known better.
The radio system thus uses four frequencies in each direction: two for Doppler
tracking, one for ranging, and one for digital communication. Up- and downlink
frequencies are slightly offset from each other to isolate the receivers at
each end from the transmitters. Each of the up/down frequency pairs has 11400
antennae in the array.

The transmitter powers of the array and the ground station are determined by
their mutual distance $d$, the ground and array antenna gains $G_{g}$ and
$G_{a}$, and the receiver sensitivities. On average, half of the 11400
antennae in the array are in view of the Earth at a given time, and their
random orientations introduce Adams (2007) another factor of $1/4$, so
$G_{a}=10\log(11400/8)=31.5\mbox{ dB}$ for all wavelengths. Signal diminuition
with distance is described by the free-space loss $L=(4\pi d/l)^{2}$, i.e.,
the reciprocal solid angle subtended by a dipole antenna at one end of the
radio link, as seen from the other. For a maximum operating distance of 1.82
AU at a frequency of 8 GHz, ${\cal L}=10\log_{10}L=-279\mbox{ dB}$. Finally,
the receiver on Earth is assumed to use a 70-meter dish antenna (gain
$G_{g}=65\mbox{ dBi}$ at 8 GHz) of the deep-space network Miller (2006). The
sum $G_{a}+G_{g}+{\cal L}$ gives a total transmission loss of -182.5 dB for
all radio frequencies (${\cal L}$ and $G_{g}$ go opposite in frequency,
$G_{a}$ is frequency-independent). The receiver system noise temperature of
the ground station is estimated Manshadi (2000) at 17 K, corresponding to a
spectral power density of -186.3 dBm/Hz. For downlink Doppler tracking, this
noise power needs to be matched within a scintillation-broadened bandwidth of
0.01 Hz. With a 20-dB noise margin and a transmission loss of -182.5 dB, the
array thus needs to output -3.8 dBm (0.4 mW) for each of the two radio
frequencies. The design should foresee a slightly higher power to allow for a
reduction in array gain due to a potential loss of radio units. For downlink
ranging at 20-MHz bandwidth and 0-dB noise margin, the transmitter pulse power
is 69 dBm (8.3 kW), to be shared among – on average – 5472 antennae. Each
antenna driver thus outputs a peak power of about 1.5 W, i.e., about as much
as a cell phone. With 3-$\mu$s-long ranging pulses transmitted once in 30 ms,
the average power for ranging is 800mW for the array, 3.4 mW per radio unit
(assuming 235 of 475 being in view of Earth and active), and 140 $\mu$W per
antenna. The lack of noise margin is compensated by repeated ranging pulses:
10000 pulses for a 20dB margin taking 300 s define the time for one ranging
measurement. Once the range is known to tens of meters, shorter PRNs can be
used for lower average power or improved noise margin. Loss of some radio
units does not disable ranging, but leads to an increase in the measurement
time.

For the uplink, a higher transmitter power is available, but the receiver
noise temperature is much higher (assumed here 1000 K , which is typical for
cell phones). Thus, the ground station has to transmit at 60 times the
downlink power, i.e., bursts of about 0.5 MW for ranging and a continuous
signal of 24 mW for Doppler tracking. The initial-phasing signal has to be
received without array gain ($G_{a}=0\mbox{ dBi}$), and it needs to have a
bandwidth given by the Doppler-shift uncertainty $\Delta\nu$ due to an initial
velocity uncertainty $\Delta v$. For $\Delta v=1m/s$ and $\nu=8$ GHz,
$\Delta\nu=27Hz$, and thus the ground station has to transmit at 1.6 kW for a
noise margin of 20 dB.

Power is generated by radiation-tolerant, thin-film, amorphous-silicon solar
cells deposited on the plastic foils of the nodes (see Fig. 4). This
technology achieves a power of 2 kW per kg of solar-cell/substrate material
United Solar Ovonic (2007) for operation at 1 A.U. distance from the sun.
Outside of the solar-cell areas, the foils have a broadband dielectric mirror
coating to reflect visible and near-infrared sunlight, but let far-infrared
escape. This reduces the temperature under the coating to well below the 321 K
in black-body radiation equilibrium at 1 A.U. from the sun, so the radio unit
can shed excess heat.

Figure 4: Left: The plastic foil of one node with 25 transponders sitting on
it. Lower right: Closeup of the foil surrounding one transponder with
microwave strip lines going out to the antennae, and a section of the plastic
foil coated with a dielectric-mirror multilayer that reflects visible and
near-infrared sunlight, but lets far infrared through. Upper right: Side view
of the same. Sunlight is reflected off the coating, but far infrared heat
radiation can escape into space. The transponder can thus radiate heat off to
the ground surrounding it.

At any given time, about half the radio units are in the sunlight and have to
power themselves and the other half through the metal wires between and within
the nodes, which are switched into rapidly varying configurations of closed
circuit loops containing power sources and sinks. At an efficiency of 10 %,
the radio transmitters will consume an average power of 50 mW (10 times (3.4
mW for ranging, 2 $\mu$W for Doppler, 1.6 mW for digital communincation), 50
mW are needed for the onboard digital electronics, and 900 mW for cooling
radio units in the sunlight, so a radio unit in the sunlight needs to harvest
1 W for itself. Radio units in the shade need 100 mW for the electronics, 100
mW power reserve for heating, and 100 mW for internal power conversion losses.
As estimated below, a sunlit radio unit then needs to harvest another 2.4 W to
power up to two units in the shade. With wires made of a high-strength
aluminum alloy with a specific resistance of $2.9\cdot 10^{-8}$ $\Omega$m
Adams (2007) and a cross section of $2.0\cdot 10^{-8}$ m2 (250 wires of 10
$\mu$m diameter, each), the wires have a resistance of 1.48 $\Omega$/m. For
example, in a 1250-meter-long current loop with two sunlit nodes powering four
in the shade (each with 25 radio units), a current of 100 mA needs to run 50 %
of the time to supply 30W at a voltage drop of 600 V (150 V per node). There
is an additional voltage drop of 185 V in the wires. Therefore, each of the
two sunlit nodes need to generate a voltage of 393 V to power the loop. At 66
% power-conversion efficiency each sunlit node has to harvest about 60 W for
those in the shade, i.e., 2.4 W per radio unit. For more detail, see ref.
Adams (2007).

To determine the payload to be delivered to the asteroid, the following masses
are estimated Adams (2007): 47.5 kg for the 475 radio, 5 kg for the plastic
foil of the 19 nodes (excluding the solar cells), 0.81 kg of solar cells
(475*3.4 W / (2 kW/kg)), 1.1 kg for 12 km of carbon/glass fiber/metal wires, 6
kg for the power-conversion electronics, master clocks and other
infrastructure in the nodes, and 9.5 kg for the containers in which the nodes
are stowed until they are deployed (Fig. 2). This brings the total array mass
to about 70 kg. Another 80 kg can be estimated for the delivery spacecraft
(structure, engines, tanks, etc.) for a mass $m_{s}=150$ kg to be delivered to
the asteroid. An optimal transfer trajectory from low-earth orbit (LEO) to
Apophis requires a total velocity change of $\Delta v=7.71$ km/s Adams (2007).
Applying the Tsiolkovsky rocket equation $m_{t}=m_{s}\exp(\Delta v/v_{p})$ for
the total mass $m_{t}$ at LEO, where $v_{p}$ is the rocket exhaust velocity
(3.07 km/s for hydrazine thrusters), the mass to launch to LEO is about 1900
kg, which is well within the capabilities of commercial satellite launchers.

##  IV Other Possible Applications

Space-based phased arrays – based on asteroids, or floating freely – could
also be used for other scientific applications, especially so if multiple ones
are used. These applications make use of the capability of a phased array to
transmit and receive simultaneously in different directions. One example is
very-long-baseline (VLBI) radio astronomy Adams (2007) with baselines the size
of Earth’s orbit. The arrays would maintain radio links with each other to
distribute a master clock signal, and would simultaneously receive
directionally from the radio source of interest. As the asteroids progress
along their orbits, different Fourier components of the source’s angular
distribution on the sky can be resolved to assemble an extremely high-
resolution image, corresponding to the spacing between the asteroid-based
receivers. Another possible application is in gravitational-wave detection at
milli-Hertz frequencies, as is currently being done with the Cassini
spacecraft Asmar et al. (2005). Because this application requires a very high
velocity accuracy, the arrays have to be based on asteroids to minimize their
susceptibility to orbital perturbations due to the solar wind, the Yarkovsky
effect, etc. With phased arrays based on several asteroids, multiple coherent
links can be maintained similtaneously in different directions between
asteroids to precisely locate a source Adams (2008). Yet another possible
application would be in the mapping of the solar wind by measuring the
scintillations in radio links between several arrays. The transmitter power
for the latter two applications may be estimated as before for the downlink in
Doppler tracking, with the bandwidth again given by scintillations. However,
the receiver temperature is now 1000K (instead of 17 K), and the receiver
antenna gain is 31.5 dBi (instead of 65 dBi). This requires an increase of the
transmitter power by a factor of $1.1\cdot 10^{5}$ over the value of 0.4 mW
estimated before. Thus, for each link the array must transmit at 45 W, and
each of the about 200 to 250 radio units in view of the respective receiver
must transmit about 200 mW. These scientific applications are not subject to
the time constraints of a mission to track an asteroid before an impact or a
keyhole passage. Therefore, lengthier but more energy-efficient transfer
trajectories Marsden and Ross (2006) are admissible, which will drastically
reduce the mass to launch to LEO, and thus the mission cost.

##  V Summary

In summary, a way is proposed to precisely track the orbit of an asteroid by
deploying a radio transponder on its surface, which establishes a two-
frequency noncoherent link for Doppler tracking at 1 $\mu$m/s accuracy,
returns broadband-modulated pulses for ranging to 5 m, and communicates
digital data representing correlations between uplink and downlink Doppler and
ranging signals. The transponder is a phased array of 11400 antennae for each
of the four frequencies used, driven by 475 radio units. These are linked with
each other by a net that serves the multiple purposes of enabling the landing
maneuver in the low gravity of the asteroid, phasing the array, and sharing
solar-generated power. Due to a high degree of redundancy, the array can
tolerate some damage, which will only lead to a slight reduction in antenna
gain. The same phased-array technology can also be used for very-long-baseline
radio astronomy or for gravitational-wave detection.

## References

  * (1) _The order of magnitude can be estimated by converting Kepler’s 3rd law to express the orbital period dependent on the velocity, and then taking the derivative._
  * Schweickart et al. (2006) R. Schweickart, C. Chapman, D. Durda, P. Hut, B. Bottke, and D. Nesvorny, White paper no. 39 at NASA NEO Workshop held in Vail, Colorado, 26-28 June, 2006 (2006), www.b612foundation.org/papers/wpdynamics.pdf.
  * Asmar et al. (2005) S. Asmar, J. Armstrong, L. Iess, and P. Tortora, Radio Science 40, 1 (2005).
  * Chesley et al. (2003) S. Chesley, S. Ostro, D. Vokrouhlicky, D. Capek, J. Giorgini, M. Nolan, J. Margot, A. Hine, L. Benner, and A. Chamberlin, Science 302, 1739 (2003).
  * Spitale (2002) J. N. Spitale, Science 296, 77 (2002).
  * Adams (2007) B. W. Adams, contribution to the Apophis mission design competition of the Planetary Society (2007), supplementary material.
  * Woo (1978) R. Woo, Astrophys. J. 219, 727 (1978).
  * Ho et al. (2002) C. Ho, M. Sue, A. Bedrossian, and R. Sniffin, Proc. GA02 1278-1287 14, 1277 (2002).
  * Armstrong (2006) J. Armstrong, Living Reviews in Relativity 9, 1 (2006).
  * Armstrong et al. (2003) J. Armstrong, L. Iess, P. Tortora, and B. Bertotti, Astrophys. J. 599, 806 (2003).
  * Howard et al. (1974) H. Howard, G. Tyler, G. Fjeldbo, A. Kliore, G. Levy, D. Brunn, R. Dickinson, R. Edelson, W. Martin, R. Postal, et al., Science 183, 1297 (1974).
  * Woo et al. (1976) R. Woo, F. Yang, K. Yip, and W. Kendall, Astrophys. J. 210, 568 (1976).
  * DeBoy et al. (2005) C. C. DeBoy, C. Haskins, D. Duven, R. Schulze, J. R. Jensen, M. Bernacik, and W. Millard, Acta Astronautica 57, 540 (2005).
  * Agilent Technologies (2005) Agilent Technologies, web document (2005), http://eesof.tm.agilent.com/docs/adsdoc2005A/gsm/gm013.html.
  * Miller (2006) R. Miller, web document (2006), www.ioag.org/ioag6_nasajpl_dsn_status.pdf.
  * Manshadi (2000) F. Manshadi, IEEE Aerospace Conf. Proc. 5, 73 (2000).
  * United Solar Ovonic (2007) L. United Solar Ovonic, web site (2007), http://www.uni-solar.com.
  * Adams (2008) B. W. Adams, manuscript in preparation (2008).
  * Marsden and Ross (2006) J. Marsden and S. Ross, Bull. Am. Amth. Soc. 43, 43 (2006).


$(m\times m)$-dimensional identity matrix as $n\rightarrow\infty$, which
follows from the assertion in Theorem 3.5 and in view of the assertion in
Theorem 3.4, the result follows. ∎

## References

  * Acemoglu et al. [2012] Daron Acemoglu, Simon Johnson, and James A. Robinson.  The colonial origins of comparative development: An empirical investigation: Reply.  _The American Economic Review_ , 102(6):3077–3110, 2012.  ISSN 00028282.
  * Angrist et al. [2006] Joshua Angrist, Victor Chernozhukov, and Iván Fernández-Val.  Quantile Regression under Misspecification, with an Application to the U.S. Wage Structure.  _Econometrica_ , 74(2):539–563, March 2006.
  * Beyhum et al. [2024] Jad Beyhum, Lorenzo Tedesco, and Ingrid Van Keilegom.  Instrumental variable quantile regression under random right censoring.  _Econom. J._ , 27(1):21–36, 2024.  ISSN 1368-4221,1368-423X.  doi: 10.1093/ectj/utad015.  URL https://doi.org/10.1093/ectj/utad015.
  * Bickel [1965] Peter J. Bickel.  On some robust estimates of location.  _The Annals of Mathematical Statistics_ , 36(3):847–858, 1965.  ISSN 00034851.
  * Blundell and Powell [2007] Richard Blundell and James L. Powell.  Censored regression quantiles with endogenous regressors.  _Journal of Econometrics_ , 141(1):65–83, 2007.  ISSN 0304-4076.  doi: https://doi.org/10.1016/j.jeconom.2007.01.016.  Semiparametric methods in econometrics.
  * Bowden Roger J. [Roger John] D. A. (Darrell A.) Bowden Roger J. (Roger John) Turkington.  _Instrumental variables / Roger J. Bowden and Darrell A. Turkington_.  Cambridge University Press, 1984.  ISBN 0521385822.
  * Camponovo [2015] Taisuke Camponovo, Lorenzo; Otsu.  Robustness of bootstrap in instrumental variable regression.  _Econometric Reviews_ , 32:352–393., 03 2015.
  * Chaudhuri [1991] Probal Chaudhuri.  Nonparametric estimates of regression quantiles and their local Bahadur representation.  _Ann. Statist._ , 19(2):760–777, 1991.  ISSN 0090-5364,2168-8966.  doi: 10.1214/aos/1176348119.  URL https://doi.org/10.1214/aos/1176348119.
  * Chen [2018] Songnian Chen.  Sequential estimation of censored quantile regression models.  _Journal of Econometrics_ , 207(1):30–52, 2018.  ISSN 0304-4076.  doi: https://doi.org/10.1016/j.jeconom.2018.06.020.
  * Chen and Wang [2020] Songnian Chen and Qian Wang.  Semiparametric estimation of a censored regression model with endogeneity.  _J. Econometrics_ , 215(1):239–256, 2020.  ISSN 0304-4076,1872-6895.  doi: 10.1016/j.jeconom.2019.08.006.  URL https://doi.org/10.1016/j.jeconom.2019.08.006.
  * Chernozhukov and Hansen [2006] Victor Chernozhukov and Christian Hansen.  Instrumental quantile regression inference for structural and treatment effect models.  _Journal of Econometrics_ , 132:491–525, 06 2006.  doi: 10.1016/j.jeconom.2005.02.009.
  * Chernozhukov and Hong [2002] Victor Chernozhukov and Han Hong.  Three-step censored quantile regression and extramarital affairs.  _Journal of the American Statistical Association_ , 97(459):872–882, 2002.  ISSN 01621459.
  * Chernozhukov et al. [2015] Victor Chernozhukov, Iván Fernández-Val, and Amanda E. Kowalski.  Quantile regression with censoring and endogeneity.  _Journal of Econometrics_ , 186(1):201–221, 2015.  ISSN 0304-4076.  doi: https://doi.org/10.1016/j.jeconom.2014.06.017.
  * Dette and Volgushev [2008] Holger Dette and Stanislav Volgushev.  Non-crossing non-parametric estimates of quantile curves.  _J. R. Stat. Soc. Ser. B Stat. Methodol._ , 70(3):609–627, 2008.  ISSN 1369-7412,1467-9868.  doi: 10.1111/j.1467-9868.2008.00651.x.  URL https://doi.org/10.1111/j.1467-9868.2008.00651.x.
  * Dette et al. [2011] Holger Dette, Jens Wagener, and Stanislav Volgushev.  Comparing conditional quantile curves.  _Scand. J. Stat._ , 38(1):63–88, 2011.  ISSN 0303-6898,1467-9469.  doi: 10.1111/j.1467-9469.2010.00718.x.  URL https://doi.org/10.1111/j.1467-9469.2010.00718.x.
  * Dhar and Shalabh [2022] Subhra Sankar Dhar and Shalabh.  GIVE statistic for goodness of fit in instrumental variables models with application to COVID data.  _Nature Scientific Reports_ , page 9472, 2022.  doi: 10.3150/22-BEJ1509.  URL https://doi.org/10.1101/2021.04.18.440376.
  * Dhar and Wu [2023] Subhra Sankar Dhar and Weichi Wu.  Comparing time varying regression quantiles under shift invariance.  _Bernoulli_ , 29(2):1527–1554, 2023.  ISSN 1350-7265,1573-9759.  doi: 10.3150/22-bej1509.  URL https://doi.org/10.3150/22-bej1509.
  * Dhar et al. [2022] Subhra Sankar Dhar, Prashant Jha, and Prabrisha Rakshit.  The trimmed mean in non-parametric regression function estimation.  _Theory Probab. Math. Statist._ , pages 133–158, 2022.  ISSN 0094-9000,1547-7363.  doi: 10.1090/tpms/1174.  URL https://doi.org/10.1090/tpms/1174.
  * Frandsen [2015] Brigham R. Frandsen.  Treatment effects with censoring and endogeneity.  _J. Amer. Statist. Assoc._ , 110(512):1745–1752, 2015.  ISSN 0162-1459,1537-274X.  doi: 10.1080/01621459.2015.1017577.  URL https://doi.org/10.1080/01621459.2015.1017577.
  * Gutenbrunner and Jureckova [1992] C. Gutenbrunner and J. Jureckova.  Regression rank scores and regression quantiles.  _The Annals of Statistics_ , 20(1):305–330, 1992.  ISSN 00905364.
  * Hampel et al. [1986] Frank R. Hampel, Elvezio M. Ronchetti, Peter J. Rousseeuw, and Werner A. Stahel.  _Robust statistics_.  Wiley Series in Probability and Mathematical Statistics: Probability and Mathematical Statistics. John Wiley & Sons, Inc., New York, 1986.  ISBN 0-471-82921-8.  The approach based on influence functions.
  * Hong and Tamer [2003] Han Hong and Elie Tamer.  Inference in censored models with endogenous regressors.  _Econometrica_ , 71(3):905–932, 2003.  ISSN 00129682, 14680262.
  * Honoré and Hu [2004] Bo E Honoré and Luojia Hu.  Estimation of cross sectional and panel data censored regression models with endogeneity.  _Journal of Econometrics_ , 122(2):293–316, 2004.  ISSN 0304-4076.  doi: https://doi.org/10.1016/j.jeconom.2003.06.001.
  * Hössjer and Karlsson [2023] Ola Hössjer and Måns Karlsson.  On the use of l-functionals in regression models.  _Open Mathematics_ , 21(1):20220597, 2023.  doi: doi:10.1515/math-2022-0597.  URL https://doi.org/10.1515/math-2022-0597.
  * Khan and Tamer [2009] Shakeeb Khan and Elie Tamer.  Inference on endogenously censored regression models using conditional moment inequalities.  _Journal of Econometrics_ , 152(2):104–119, 2009.  ISSN 0304-4076.  doi: https://doi.org/10.1016/j.jeconom.2009.01.006.  Nonparametric and Robust Methods in Econometrics.
  * Kim et al. [2011] Daniel Kim, Christopher Baum, Michael Ganz, S Subramanian, and Ichiro Kawachi.  The contextual effects of social capital on health: A cross-national instrumental variable analysis.  _Social science and medicine (1982)_ , 73:1689–97, 10 2011.  doi: 10.1016/j.socscimed.2011.09.019.
  * Knight [1998] Keith Knight.  Limiting distributions for $L_{1}$ regression estimators under general conditions.  _Ann. Statist._ , 26(2):755–770, 1998.  ISSN 0090-5364,2168-8966.  doi: 10.1214/aos/1028144858.  URL https://doi.org/10.1214/aos/1028144858.
  * Koenker and Bassett [1978] Roger Koenker and Gilbert Bassett, Jr.  Regression quantiles.  _Econometrica_ , 46(1):33–50, 1978.  ISSN 0012-9682,1468-0262.  doi: 10.2307/1913643.  URL https://doi.org/10.2307/1913643.
  * Koenker and Bassett [1982] Roger Koenker and Gilbert Bassett, Jr.  Robust tests for heteroscedasticity based on regression quantiles.  _Econometrica_ , 50(1):43–61, 1982.  ISSN 0012-9682,1468-0262.  doi: 10.2307/1912528.  URL https://doi.org/10.2307/1912528.
  * Koenker and Portnoy [1987] Roger Koenker and Stephen Portnoy.  L-estimation for linear models.  _Journal of the American Statistical Association_ , 82(399):851–857, 1987.  ISSN 01621459.
  * Lee [1992] Myoung-Jae Lee.  Winsorized mean estimator for censored regression.  _Econometric Theory_ , 8(3):368–382, 1992.  ISSN 02664666, 14694360.
  * Lee [2021] Myoung-jae Lee.  Instrument residual estimator for any response variable with endogenous binary treatment.  _J. R. Stat. Soc. Ser. B. Stat. Methodol._ , 83(3):612–635, 2021.  ISSN 1369-7412,1467-9868.  doi: 10.1111/rssb.12442.  URL https://doi.org/10.1111/rssb.12442.
  * Lu and Li [2021] Jiarui Lu and Hongzhe Li.  Hypothesis testing in high-dimensional instrumental variables regression with an application to genomics data.  _Statist. Sinica_ , 32:613–633, 2021.  ISSN 1017-0405,1996-8507.  doi: 10.1007/s11045-020-00752-x.  URL https://doi.org/10.1007/s11045-020-00752-x.
  * Newey [1987] Whitney K. Newey.  Efficient estimation of limited dependent variable models with endogenous explanatory variables.  _Journal of Econometrics_ , 36(3):231–250, 1987.  ISSN 0304-4076.  doi: https://doi.org/10.1016/0304-4076(87)90001-7.
  * Powell [1984] James L Powell.  Least absolute deviations estimation for the censored regression model.  _Journal of Econometrics_ , 25(3):303–325, 1984.  ISSN 0304-4076.  doi: https://doi.org/10.1016/0304-4076(84)90004-6.
  * Powell [1986a] James L. Powell.  Censored regression quantiles.  _Journal of Econometrics_ , 32(1):143–155, 1986a.  ISSN 0304-4076.  doi: https://doi.org/10.1016/0304-4076(86)90016-3.
  * Powell [1986b] James L. Powell.  Symmetrically trimmed least squares estimation for tobit models.  _Econometrica_ , 54(6):1435–1460, 1986b.  ISSN 00129682, 14680262.
  * Rivers and Vuong [1988] Douglas Rivers and Quang H. Vuong.  Limited information estimators and exogeneity tests for simultaneous probit models.  _Journal of Econometrics_ , 39(3):347–366, 1988.  ISSN 0304-4076.  doi: https://doi.org/10.1016/0304-4076(88)90063-2.
  * Serfling [1980] Robert Serfling.  _Approximation Theorems of Mathematical Statistics_.  Wiley, 1980.  URL https://api.semanticscholar.org/CorpusID:122686903.
  * Shalabh and Dhar [2023] Shalabh and Subhra Sankar Dhar.  Testing the goodness of fit in instrumental variables models.  In _G Families of Probability Distributions_ , pages 330–343. Taylor & Francis, 2023.
  * Shukla et al. [2023] Swati Shukla, Subhra Sankar Dhar, and Shalabh.  M-estimation in censored regression model using instrumental variables under endogeneity, 2023.
  * Smith and Blundell [1986] Richard J. Smith and Richard W. Blundell.  An exogeneity test for a simultaneous equation tobit model with an application to labor supply.  _Econometrica_ , 54(3):679–685, 1986.  ISSN 00129682, 14680262.
  * Taylor and Otsu [2019] Luke Taylor and Taisuke Otsu.  Estimation of nonseparable models with censored dependent variables and endogenous regressors.  _Econometric Rev._ , 38(1):4–24, 2019.  ISSN 0747-4938,1532-4168.  doi: 10.1080/07474938.2016.1235310.  URL https://doi.org/10.1080/07474938.2016.1235310.
  * van der Vaart and Wellner [1996a] Aad W. van der Vaart and Jon A. Wellner.  _Weak convergence and empirical processes_.  Springer Series in Statistics. Springer-Verlag, New York, 1996a.  ISBN 0-387-94640-3.  doi: 10.1007/978-1-4757-2545-2.  URL https://doi.org/10.1007/978-1-4757-2545-2.  With applications to statistics.
  * van der Vaart and Wellner [1996b] Aad W. van der Vaart and Jon A. Wellner.  _Independence Empirical Processes_ , pages 367–371.  Springer New York, New York, NY, 1996b.  ISBN 978-1-4757-2545-2.  doi: 10.1007/978-1-4757-2545-2˙34.  URL https://doi.org/10.1007/978-1-4757-2545-2_34.
  * Vansteelandt and Didelez [2018] Stijn Vansteelandt and Vanessa Didelez.  Improving the robustness and efficiency of covariate‐adjusted linear instrumental variable estimators.  _Scandinavian Journal of Statistics_ , 45:941 – 961, 2018.
  * Vapnik and Červonenkis [1971] V. N. Vapnik and A. Ja. Červonenkis.  The uniform convergence of frequencies of the appearance of events to their probabilities.  _Teor. Verojatnost. i Primenen._ , 16:264–279, 1971.  ISSN 0040-361x.
  * Wang and Zhang [2024] Wenjie Wang and Yichong Zhang.  Wild bootstrap inference for instrumental variables regressions with weak and few clusters.  _J. Econometrics_ , 241(1):Paper No. 105727, 2024.  ISSN 0304-4076,1872-6895.  doi: 10.1016/j.jeconom.2024.105727.  URL https://doi.org/10.1016/j.jeconom.2024.105727.
  * White [1980] Halbert White.  Nonlinear regression on cross-section data.  _Econometrica_ , 48(3):721–746, 1980.  ISSN 00129682, 14680262.
  * Zhou and Shao [2013] Zhou Zhou and Xiaofeng Shao.  Inference for linear models with dependent errors.  _Journal of the Royal Statistical Society. Series B (Statistical Methodology)_ , 75(2):323–343, 2013.  ISSN 13697412, 14679868.