# Learned 1-D passive scalar advection to accelerate chemical transport
modeling: a case study with GEOS-FP horizontal wind fields

###### Abstract

We developed and applied a machine-learned discretization for one-dimensional
(1-D) horizontal passive scalar advection, which is an operator component
common to all chemical transport models (CTMs). Our learned advection scheme
resembles a second-order accuracy, three-stencil numerical solver, but differs
from a traditional solver in that coefficients for each equation term are
output by a neural network rather than being theoretically-derived constants.
We downsampled higher-resolution simulation results—resulting in up to
16$\times$ larger grid size and 64$\times$ larger timestep—and trained our
neural network-based scheme to match the downsampled integration data. In this
way, we created an operator that is low-resolution (in time or space) but can
reproduce the behavior of a high-resolution traditional solver. Our model
shows high fidelity in reproducing its training dataset (a single 10-day 1-D
simulation) and is similarly accurate in simulations with unseen initial
conditions, wind fields, and grid spacing. In many cases, our learned solver
is more accurate than a low-resolution version of the reference solver, but
the low-resolution reference solver achieves greater computational speedup
(500$\times$ acceleration) over the high-resolution simulation than the
learned solver is able to (18$\times$ acceleration). Surprisingly, our learned
1-D scheme—when combined with a splitting technique—can be used to predict 2-D
advection, and is in some cases more stable and accurate than the low-
resolution reference solver in 2-D. Overall, our results suggest that learned
advection operators may offer a higher-accuracy method for accelerating CTM
simulations as compared to simply running a traditional integrator at low
resolution.

Chemical transport modeling (CTM) is an essential tool to study air pollution.
CTM simulation takes a long computing time. Modeling pollutant transport
(advection) is the second most computationally intensive part of the model.
Decreasing the resolution reduces the advection computing time, but decreases
accuracy. We employed machine learning to reduce the resolution of advection
while keeping the accuracy. We verified the robustness of our solver with
several generalization testing scenarios. In our 2-D simulation, our solver
showed a 100-time faster simulation with fair accuracy. Integrating our
approach to existing CTM will enable much more simulations than now to study
air pollution and to provide solutions.

##  1 Introduction

Atmospheric chemical transport models (CTMs) are a key research tool to study
air pollution and to predict outcomes of pollution mitigation efforts.
Currently available CTMs (e.g., GEOS-Chem (Bey et al., 2001), CMAQ (EPA,
2022), and WRF-Chem (Grell et al., 2005)) numerically represent
physicochemical processes including emission, chemical reactions, transport,
and deposition discretized across many grid boxes. This numerical
discretization engenders a dilemma between accuracy and computational cost:
higher-resolution simulations are more accurate but can be computationally
intractable for some use cases. Parallel and distributed computation can
reduce the overall execution time of simulations but add cost and complexity.
For example, Zhuang et al. (2020) found that an $8\times$ increase in
distributed computing power resulted in only a $4\times$ increase in speed for
GEOS-Chem simulations.

Recent efforts have explored the use of machine-learned (ML) surrogate
modeling to address computational costs in CTMs. Because the chemistry solver
is the most expensive part of the atmospheric chemistry model, most ML
research related to the acceleration of CTMs has focused on the chemistry
operator (Keller et al., 2017; Keller and Evans, 2019; Shen et al., 2022; Kelp
et al., 2020, 2022; Huang and Seinfeld, 2022). The advective transport
operator is the second-most expensive module of the air quality model (Eastham
et al., 2018), but we are only aware of one existing study of machine-learned
advection (Zhuang et al., 2021). Sturm et al. (2023) accelerated the advection
using by compressing the number of species, rather than accelerating the
advection solver itself. However, recent progress in machine-learned
computational fluid dynamics (CFD) and turbulence simulation (Kochkov et al.,
2021; Stachenfeld et al., 2021; Brunton et al., 2020; Vinuesa and Brunton,
2022), has laid a solid foundation for the acceleration of the transport
operator in the context of air quality modeling.

In this study, we develop a machine-learned 1-D horizontal advection solver in
the context of an air quality model. We build on the work of Zhuang et al.
(2021), who trained a model using synthetic velocity fields within the
rectangular computational domain with a fixed grid size. Here, as a step
toward full-scale deployment in an ML advection operator in a CTM, we instead
use real-world historical wind velocity data and grid cells that change size
with changing latitude. Furthermore, we extensively explore the robustness of
our approach by limiting the training data used and by testing its performance
under a variety of different conditions. Altogether, the results herein
characterize the promise of machine-learned advection operators for use in
CTMs and outline the remaining challenges to their full-scale adoption.

The remaining parts of this paper are structured as follows: Section 2
describes our methodology, including the numerical advection scheme used as a
reference solver, the structure of the learned scheme which we train to match
the results of our reference solver, and numerical experiments including
generalization testing and 2-D application. Section 3 presents the results of
the experiments introduced in the methodology section. Section 4 discusses
strengths and limitations of our approach, and future research needs. Section
5 concludes this paper.

##  2 Methodology

###  2.1 Numerical advection

####  2.1.1 Numerical scheme and dataset generation

We implemented the so-called L94 advection scheme (Lin et al., 1994) using the
Julia scientific computing language (Bezanson et al., 2012). The L94 advection
scheme is a second-order-accurate van Leer-type advection scheme and has been
used in a well-established chemical transport model, GEOS-Chem (Bey et al.,
2001), when coupled with a nondirectional multidimensional splitting (Lin and
Rood, 1996). An implementation in Julia—a differentiable programming language
(Innes et al., 2019)—allows us to train machine-learned models by
backpropagating error gradients through multiple model time steps. Following
Lin et al. (1994), the equations for our 1-D numerical advection reference
model are as follows:

Suppose we have a scalar field discretized across a 1-D grid, with
concentration $\Phi_{i}^{n}$ in the ith grid point at the $n^{\rm th}$
timestep. Then, the cell average value ($\Phi_{i+1/2}^{n}$) can be computed as
Equation 1 because the L94 scheme assumes the piecewise linear distribution of
the scalar field inside a cell:

$\Phi_{i+1/2}^{n}=(\Phi_{i}^{n}+\Phi_{i+1}^{n})/2.$ (1)

The ultimate purpose of this flux form operator is to calculate the next time
cell average as in Equation 2:

$\Phi_{i+1/2}^{n+1}=\Phi_{i+1/2}^{n}-\frac{\Delta t}{\Delta
x}[f_{(i+1)}-f_{(i)}],$ (2)

where the flux ($f_{(i)}$) through the cell boundary can be expressed by
Equation 3a for $U_{i}\geq$ 0 or Equation 3b for $U_{i}<0$, where $U_{i}$ is
the velocity at the ith edge of the grid cell:

$\displaystyle
f_{(i)}=U_{i}[\Phi_{i-1/2}^{n}+\frac{\Delta\Phi_{i-1/2}^{n}}{2}(1-C_{i}^{-})],$
(3a) $\displaystyle
f_{(i)}=U_{i}[\Phi_{i+1/2}^{n}+\frac{\Delta\Phi_{i+1/2}^{n}}{2}(1+C_{i}^{+})].$
(3b)

Both $C_{i}^{-}$ and $C_{i}^{+}$ are Courant–Friedrichs–Lewy (CFL) numbers
(Courant et al., 1928), which can be calculated using Equation 4a for
$U_{i}<0$ and Equation 4b for $U_{i}\geq$ 0, respectively. Here $\Delta
x_{i+-1/2}$ and $\Delta x_{i+1/2}$ are the grid spacing on the left and right
sides of the $i^{\rm th}$ grid cell edge, and $\Delta t$ is the timestep.

$\displaystyle C_{i}^{-}=\frac{U_{i}\Delta t}{\Delta x_{i-1/2}},$ (4a)
$\displaystyle C_{i}^{+}=\frac{U_{i}\Delta t}{\Delta x_{i+1/2}}.$ (4b)

The key feature of this reference scheme is the method for deriving the
derivative of the cell-averaged scalar, which is implemented using a
monotonicity constraint (denoted “mono5” following Lin et al. (1994)) to
ensure numerical stability as shown in Equation 5:

$[\Delta\Phi_{i+1/2}]_{\textnormal{mono5}}=\textnormal{sign}([\Delta\Phi_{i+1/2}]_{\textnormal{avg}})\times\textnormal{min}[|[\Delta\Phi_{i+1/2}]_{\textnormal{avg}}|,2(\Phi_{i+1/2}-\Phi_{i+1/2}^{\textnormal{min}}),2(\Phi_{i+1/2}^{\textnormal{max}}-\Phi_{i+1/2})].$
(5)

In Equation 5, $\Phi_{i+1/2}^{\textnormal{max}}$ and
$\Phi_{i+1/2}^{\textnormal{min}}$ represent the local maximum and minimum,
respectively, of $\Phi$ in the three-point stencil centered at $i+1/2$. Then,
$[\Delta\Phi_{i+1/2}]_{\textnormal{avg}}$ is the average of the local spatial
difference in $\Phi$ as in Equation 6:

$[\Delta\Phi_{i+1/2}]_{\textnormal{avg}}=\frac{(\delta\Phi_{i}+\delta\Phi_{i+1})}{2},$
(6)

where

$\delta\Phi_{i}=\Phi_{i+1/2}^{n}-\Phi_{i-1/2}^{n}.$ (7)

Our 1-D simulation domain is a straight line across the GEOS-FP
$0.25\text{\,}\mathrm{\SIUnitSymbolDegree}$ $\times$
$0.3125\text{\,}\mathrm{\SIUnitSymbolDegree}$ North America nested grid (NASA,
2022) which covers
$130\text{\,}\mathrm{\SIUnitSymbolDegree}$W–$60\text{\,}\mathrm{\SIUnitSymbolDegree}$W
and
$9.75\text{\,}\mathrm{\SIUnitSymbolDegree}$N–$60\text{\,}\mathrm{\SIUnitSymbolDegree}$N.
We used the component of the ground-level wind field aligned with the domain
line direction, using velocity data from
$0.25\text{\,}\mathrm{\SIUnitSymbolDegree}$ $\times$
$0.3125\text{\,}\mathrm{\SIUnitSymbolDegree}$ GEOS-FP wind fields (NASA,
2022). For example, if the line lay in the east-west direction, we would use
the east-west component of the given wind field.

We generated a training dataset by simulating advection through the East-West
horizontal line on $39.00\text{\,}\mathrm{\SIUnitSymbolDegree}$N with
$0.3125\text{\,}\mathrm{\SIUnitSymbolDegree}$ spatial resolution, using a
300-second timestep with the reference solver described above. The thick blue
line in Figure 1 shows the training domain. The simulation period is the first
10 days of January 2019. We fed a square shape initial condition which has a
100 ppb mixing ratio on the central one-third of the domain and zero on the
rest of the domain. We chose 100 ppb for the initial mixing ratio of the
passive tracer as it is a typical mixing ratio for atmospheric pollutants and
the default initial mixing ratio of the passive scalar in the GEOS-Chem
transport tracer simulation. The square-shaped initial condition results in a
challenging training dataset because the integrated wave will have a stiff
gradient at the beginning of the simulation, which becomes smoother as the
simulation progresses. We specified zero-gradient spatial boundary conditions.
This results in simulations with 224 grid cells and 2880 integration
timesteps.

Figure 1: Visualization of the 1-D training and testing domain.

####  2.1.2 Spatio-temporal coarsening

We down-sampled the scalar outputs of the advection simulation with factors of
1$\times$, 2$\times$, 4$\times$, 8$\times$ and 16$\times$ in space and
1$\times$, 2$\times$, 4$\times$, 8$\times$, 16$\times$, 32$\times$, and
64$\times$ in time. We applied both spatial and temporal down-sampling
together, resulting in 35 different scenarios of data coarsening. We down-
sampled the output of the reference simulation by averaging as shown in Figure
11. At maximum coarsening in both space and time, the data dimensionality was
reduced from 224 cells $\times$ 2,880 steps to 14 cells $\times$ 45 steps. We
then trained machine-learning models to reproduce the down-sampled reference
simulation.

###  2.2 Development of the learned advection scheme

####  2.2.1 Model structure

As illustrated in Figure 2, we designed a physics-informed machine learning
solver that is fed with scalar and velocity fields and returns the surrogate
coefficients for numerical integration of a 3-stencil second-order-accurate
scheme. We represented the temporal gradient using six terms describing
3-stencil first- and second-order derivatives as shown in Equation 8:

$\Phi_{i}^{n+1}=\Phi_{i}^{n}+k_{1}\frac{\Delta t}{\Delta
x}[a_{1},a_{2},a_{3}]\cdot[\Phi_{i-1}^{n},\Phi_{i}^{n},\Phi_{i+1}^{n}]+k_{2}\left(\frac{\Delta
t}{\Delta
x}\right)^{2}[b_{1},b_{2},b_{3}]\cdot[\Phi_{i-1}^{n},\Phi_{i}^{n},\Phi_{i+1}^{n}],$
(8)

where $k_{1}$ and $k_{2}$ are constant-gradient scaling factors (to help with
training) and $a_{i}$ and $b_{i}$ are the surrogate coefficients generated by
the learned solver. This formulation gives our learned scheme the flexibility
to use different values of $\Delta$x, which can represent the variable
relationship between grid cell size in meters and degrees at different
latitudes. Overall, Equation 8 results in a machine-learned model with a
strong inductive bias based on the mathematical framework described in Section
2.1.1.

Figure 2: Design of the neural net-based models to emulate the L94 advection
solver. $\Phi^{n}$ is the scalar field and $U^{n}$ is the velocity field.

We used a 3-layer convolutional neural network (CNN) with 32 filters per layer
to fit the surrogate coefficients. Zhuang et al. (2021) used a 4-layer CNN
with 32 filters per layer to emulate 1-D advection and a 10-layer CNN with 32
filters per layer for 2-D. We used 32 filters per layer as they did and
reduced the number of layers from 4 to 3 to reduce computational cost. We used
the Flux.jl (Innes et al., 2018) machine learning software library to
implement the neural network. We trained a separate neural network for each
combination of temporal and spatial coarsening. We used a default kernel size
of 3, corresponding to a 3-stencil scheme, but for simulations where the
maximum CFL number is larger than 1, we increased the CNN kernel size (and
stencil size) to 5, 9, 17, or 31 to provide a suitable spatiotemporal
information horizon. For example, if the maximum CFL number of the downsampled
velocity field is between 1 and 2, we set the kernel size to be five (i.e.,
two cells on the left side, one in the center, and two cells on the right
side) to allow our algorithm to access information about all scalar values
which may pass through the grid cell of interest during a given timestep. To
facilitate training our neural network across the stiff gradients that
characterize this system, we employed the Gaussian error linear unit (GeLU;
Hendrycks and Gimpel (2016)) activation function as suggested by Kim et al.
(2021). For the final activation layer, we used hyperbolic tangent to limit
the rate of potential error growth. To normalize the temporal gradient for
each term, we multiplied scaling factors on each term ($k_{1}$ and $k_{2}$ in
Equation 8), with magnitudes similar to the inverse of $\Delta t/\Delta x$ and
$(\Delta t/\Delta x)^{2}$. In this study, we used 100 for $k_{1}$, and 10,000
for $k_{2}$, because $\Delta t=300$ seconds and $\Delta x\approx 30,000$ km in
the baseline resolution.

####  2.2.2 Neural net training

We trained a separate neural network for each coarse-graining resolution
listed in Section 2.1.2. We trained our neural network to minimize the
absolute error of a 10-step prediction following Zhuang et al. (2021)’s
suggestion and to reflect the temporal variation of the advected scalar field.
(This means that each gradient descent step involves running the model forward
for 10 steps with the current neural network weights, calculating the gradient
of the average error among the 10 steps with respect to the weights, and then
updating the weights in the opposite direction of the gradient.) In our
experiments, the use of a single-timestep loss function resulted in
instability or mode collapse. Input feature standardization was performed by
multiplying the scalar array by $10^{7}$ and dividing the velocity array by 15
to ensure both scalar and velocity would be bound to a scale of $10^{0}$ for
compatibility with the default parameters of the optimizer that we use. We
introduced uniformly-distributed noise with a maximum value of $3\times
10^{-4}$ times the magnitude of the initial condition into the scalar array to
let the model learn the physical patterns behind the noise and prevent error
accumulation, as shown by Stachenfeld et al. (2021). We used the Adam
optimizer (Kingma and Ba, 2014) with the default parameters in Flux.jl. We
modified the learning rate ($\rho$) to decay with the number of training
epochs as in Equation 9:

$\rho=\frac{2\times 10^{-3}}{1+\rm epoch}.$ (9)

The purpose of the decaying learning rate is to encourage model convergence.
We trained the model for 100 epochs by default and saved the model parameters
every epoch. When the local minimum of training loss was observed within 100
epochs, we chose the training epoch as the best epoch. When the training loss
was still decreasing after 100 epochs, we continued training until the loss
stopped decreasing. The training epoch with the best performance for each
resolution is summarized in Table 1.

After training the model, we evaluated the correctness of our learned schemes
in emulating the reference solver. To do this, we fed our learned models with
the same initial conditions and integrated them with the same velocity fields
as were used for the reference simulation. We evaluated the performance of the
learned scheme using three different statistics: mean absolute error, root
mean square error, and R2, all averaged throughout the entire simulation
rather than for a single timestep. Since error tends to accumulate with time
integration, the scalar fields that the scheme faces in each step are never
exactly the same as the scalar fields fed in the training procedure.
Therefore, even though we use the same initial conditions and velocity fields
used in the training dataset, the recurrent input data in each model
simulation is somewhat distinct from the reference simulation data, so we are
not testing on the training data here, but further generalization tests are
described below.

We evaluated the computational acceleration by the learned solver against the
finest reference solver. The computational time to integrate the same duration
(10 days) by each solver was measured using a single CPU (dual 6248 Cascade
Lake CPU) within an HPE Apollo 6500 node installed in the University of
Illinois Campus Cluster. We used the $@benchmark$ macro of BenchmarkTools.jl
(JuliaCI, 2022) to collect the statistics in time evaluation. We calculated
the speedup as shown in Equation 10:

${\rm speedup}=\frac{{\rm mean}({\rm time_{learned}})}{{\rm mean}({\rm
time_{baseline}})},$ (10)

where ${\rm time_{learned}}$ is the computational time spent by the learned
solver and ${\rm time_{baseline}}$ is the computational time spent by the
reference solver in the finest resolution. When the learned solver is faster
than the baseline, the speedup is larger than 1.

###  2.3 Generalization ability testing

We evaluated the generalization ability of learned schemes against a longer
time span (3$\times$ longer than the training dataset), different initial
conditions (Dirac-delta shape and Gaussian distribution shape), and different
wind fields. The tests on the different wind fields can be further categorized
by seasons (the first 10 days in April, July, and October to represent Spring,
Summer, and Fall respectively), latitudes
($29.5\text{\,}\mathrm{\SIUnitSymbolDegree}$N and
$45\text{\,}\mathrm{\SIUnitSymbolDegree}$N), and domain direction
(longitudinal line at $76.875\text{\,}\mathrm{\SIUnitSymbolDegree}$W using the
north-south wind velocity component). The pale blue lines in Figure 1 show the
spatial domains of the generalization tests. For the longer-term stability
test, we simulated scalar advection from January 1st to 30th. For other tests,
the integration time span is 10 days as in the training data.

Each test has a different purpose. The longer-duration simulations are used to
evaluate the degradation of predictive performance over longer time horizons.
Changing the shape of the initial conditions tests whether the learned schemes
were over-fitted to the initial condition shape used in the training dataset.
The Dirac-delta initial condition test investigates the performance of the
learned scheme with extremely large spatial gradients. The purpose of using
different velocity fields is to assess the robustness of the learned scheme on
the different wind patterns. Finally, the tests with changes in latitude and
domain direction evaluate the learned schemes’ ability to work with different
grid sizes. For context, the grid size of the training set ($\Delta$x,
$0.3125\text{\,}\mathrm{\SIUnitSymbolDegree}$ in
$39\text{\,}\mathrm{\SIUnitSymbolDegree}$N) is approximately 27.0 km, while
$\Delta$x in $29.5\text{\,}\mathrm{\SIUnitSymbolDegree}$N and
$45\text{\,}\mathrm{\SIUnitSymbolDegree}$N are 24.6 and 30.3 km, respectively.
The longitudinal grid size,
$\Delta$y—$0.25\text{\,}\mathrm{\SIUnitSymbolDegree}$—is approximately 27.8 km
at all longitudes.

###  2.4 Multidimensional splitting

We implemented 2-D advection in the same way for both the reference solver and
the learned scheme developed in this study. In both cases, 1-D advection is
converted to 2-D using non-directional splitting as described in Lin and Rood
(1996) and reproduced in Equation 11:

$Q^{n+1}=Q^{n}+F(Q^{n}+\frac{1}{2}G(Q^{n}))+G(Q^{n}+\frac{1}{2}F(Q^{n}))$ (11)

where $Q^{n}$ is the scalar value in a specific point on the 2-D domain at
time $n$, and $F$ and $G$ are the x-axis and y-axis discretization operators,
respectively. This means that for the learned solver, we used models to make
2-D predictions that were only trained to make 1-D predictions (at the same
spatial and temporal resolution). All simulations were performed over the
GEOS-FP $0.25\text{\,}\mathrm{\SIUnitSymbolDegree}$ $\times$
$0.3125\text{\,}\mathrm{\SIUnitSymbolDegree}$ North America nested grid
($130\text{\,}\mathrm{\SIUnitSymbolDegree}$W–$60\text{\,}\mathrm{\SIUnitSymbolDegree}$W,
$9.75\text{\,}\mathrm{\SIUnitSymbolDegree}$N–$60\text{\,}\mathrm{\SIUnitSymbolDegree}$N).
The simulation period is the first 10 days in January 2019. The baseline
resolution is $0.25\text{\,}\mathrm{\SIUnitSymbolDegree}$ $\times$
$0.3125\text{\,}\mathrm{\SIUnitSymbolDegree}$ in space, and 5 minutes in time.
Like the 1-D study, coarse-graining was performed by doubling the spatial grid
size and the timestep several times from
$0.25\text{\,}\mathrm{\SIUnitSymbolDegree}$ $\times$
$0.3125\text{\,}\mathrm{\SIUnitSymbolDegree}$ $\times$ 5 minutes and up to
$4.0\text{\,}\mathrm{\SIUnitSymbolDegree}$ $\times$
$5.0\text{\,}\mathrm{\SIUnitSymbolDegree}$ $\times$ 5 hours 20 minutes.

##  3 Results

###  3.1 Model performance in emulating the training dataset

Figure 3 summarizes trained model performance in integrating the training
dataset with different levels of spatial and temporal coarsening. There are
three cases where the learned solver returned unstable outputs, which are
discussed in the final paragraph of this section. Except for these unstable
cases, the overall errors are fairly small. MAEs and RMSEs range from
2.47–24.8 ppb and 3.94–45.38 ppb, respectively (Figure 3A and B). MAE and RMSE
values are mostly less than 10 ppb. $R^{2}$ values are mostly higher than 0.87
with a worst case of 0.70 (Figure 3C). The maximum computational speedup
compared to the high-resolution reference model is 18.0$\times$ at the
coarsest spatial and temporal resolution. However, under many conditions, our
ML model is slower than the reference model owing to the computational
intensity of the neural network. We find that temporal coarsening reduces
computation time proportionately to the coarsening factor (e.g., 4$\Delta t$
is $2\times$ faster than 2$\Delta t$). Spatial coarsening also reduces
computation time, but sub-proportionally to the coarsening factor (e.g.,
4$\Delta x$ is $<2\times$ faster than 2$\Delta x$) because operations that are
sequential in space can leverage the single instruction multiple data (SIMD)
compiler optimizations whereas operations that are sequential in time cannot.
The statistics of the computational time by the learned solvers and the
reference solvers are given in Tables 2 and 3.

Figure 3: Performance of the learned solver in emulating the training dataset
at different coarse-graining resolutions as compared to the highest-resolution
reference solver simulation. Conditions where integration failed are marked
“Unstable”.

Computational speedup values reported in Figure 3 and elsewhere are computed
against our Julia implementation of the reference model described in Section
2.1.1. As a sense-check to demonstrate that any ML model speedups that we
observed are not owing to an overly slow implementation of the reference
model, we compare the execution time of our reference model to that of a GEOS-
Chem “transport tracer” configuration, which only computes advection. We
turned off other modules including chemistry, depositions, radiation,
convection, boundary layer mixing, etc. On a single CPU in the Illinois Campus
Cluster hardware described above, the GEOS-Chem “transport tracer”
configuration in the nested $0.25\text{\,}\mathrm{\SIUnitSymbolDegree}$
$\times$ $0.3125\text{\,}\mathrm{\SIUnitSymbolDegree}$ North America domain
with $4.0\text{\,}\mathrm{\SIUnitSymbolDegree}$ $\times$
$5.0\text{\,}\mathrm{\SIUnitSymbolDegree}$ global simulation for providing
boundary conditions takes 53,677 seconds whereas the global
$4.0\text{\,}\mathrm{\SIUnitSymbolDegree}$ $\times$
$5.0\text{\,}\mathrm{\SIUnitSymbolDegree}$ simulation takes 613 seconds for 10
days simulation. Therefore, the computational time for a 10-day simulation in
nested North America is to be 53,065 seconds, or 1.89 microseconds per grid
cell per timestep using a single CPU core, whereas our reference solver takes
0.06 microseconds per grid cell per timestep using a single CPU core. Although
this comparison does not show that our reference solver is faster than GEOS-
Chem—a “transport tracer” simulation does in fact not only do advection so it
is not a fair comparison—but it does suggest that our reference algorithm is
not unreasonably slow.

To evaluate the benefit of our learned advection model, we compare its
performance with simulations using the reference method run at spatial and
temporal resolutions matching those of the learned models. Figure 4 summarizes
the performance of the lower-resolution reference solver. Interestingly, the
coarsening with the reference solver appears to have a preferential regime
that lies on the diagonal elements in the performance matrix (e.g., 2$\Delta
x$-2$\Delta t$, 4$\Delta x$-4$\Delta t$, etc). Previous research has suggested
that shorter timesteps unconditionally result in higher model accuracy (Philip
et al., 2016), but results here may provide contradictory evidence. Our
results instead suggest that matching temporal and spatial coarsening may in
some cases enhance the performance of the model as well as the computational
speed, possibly because temporal coarsening averages out high-frequency
features, which could trigger numerical diffusion in low spatial resolution.
Also, the native timestep in the GEOS-FP wind field is 60 minutes, thus,
temporal coarsening from 5 minutes to a longer step can be beneficial by
reducing unnecessary computation and minimizing possible numerical error
accumulation. Regardless, additional investigation is required to fully
understand the phenomenon we observe here.

Figure 4: Performance of the reference solver in coarse-graining the fine
resolution training dataset in different resolutions as compared to the
highest-resolution reference solver simulation. Conditions where integration
failed are marked “Unstable”.

When both solvers are run at the same resolution, the computational speedup by
the reference solver is higher than that by the learned solver, owing to the
computational complexity of the neural network (Figures 3D and 4D). For
resolutions where we use a 3-layer CNN, the low-resolution reference solver
was 18.7–28.5 times faster than the learned solver. For the (1$\Delta x$,
4$\Delta t$) and (1$\Delta x$, 32$\Delta t$) resolutions, we add an additional
CNN layer to prevent instability, resulting in learned solvers that are
$\sim$1,000$\times$ slower than the low-resolution reference solvers. Given
the fact that the reference solver with certain coarsening factors such as
(8$\Delta x$, 64$\Delta t$) could achieve much faster speedup than the maximal
speedup from the learned solver as well as favorable accuracy, we cannot say
that the learned solver in this study is strictly superior to the traditional
solver. Instead, we would like to argue that the learned solver presented here
could be a complementary option in cases where the traditional solver does not
work well. For example, the reference solver does not work well in cases with
large grid cells combined with short timesteps (Figure 4)—such as nested
simulations with coarse global domains but small timesteps to match the
higher-resolution inner domains—but the learned solver can still operate with
high accuracy in this regime. GEOS-Chem Classic operates the transport with 10
min by default in the global simulation and 5 min by default in the nested
simulation. Using our learned solver, we can increase the time interval to
even more than 10 minutes.

The learned solvers with the highest accuracy predicting the training dataset
(8$\Delta x$, 16$\Delta t$) and largest acceleration (16$\Delta x$, 64$\Delta
t$) are shown in Figure 5, which shows that the spatial patterns with time
evolution are well-described by the learned schemes in both resolutions.
Accuracy is not perfect, and deviations are more noticeable in (16$\Delta x$,
64$\Delta t$), but $R^{2}$ values combined across all timesteps are high in
both cases (0.99 and 0.90, respectively).

Figure 5: Time series representation of the learned emulation of the training
dataset. The top row is 8$\Delta x$, 16$\Delta t$ and the bottom row is
16$\Delta x$, 64$\Delta t$. The initial shapes are slightly different from the
square initial condition because of coarse-graining.

We are unable to successfully train the model with the coarsening factors of
(1$\Delta x$, 8$\Delta t$), (1$\Delta x$, 16$\Delta t$), and (1$\Delta x$,
64$\Delta t$). In these cases, the learned model simulations become
numerically unstable. One possible explanation is that those regimes have
large CFL numbers due to large temporal coarsening factor relative to that in
space. (The low-resolution reference model also become unstable under similar
conditions.) The maximum CFL numbers in (1$\Delta x$, 8$\Delta t$), (1$\Delta
x$, 16$\Delta t$), and (1$\Delta x$, 64$\Delta t$) are 1.82, 3.63, 7.26, and
14.52, respectively. The maximum CFL numbers in all the tested coarsening
cases are summarized in Table 4. We use larger stencils and larger neural
networks in such cases to account for the increased spatio-temporal transport
horizon; this solved the problem for (1$\Delta x$, 4$\Delta t$) and (1$\Delta
x$, 32$\Delta t$) but not for the cases above. Resolving this instability
within a computationally efficient framework is a topic for future research.

###  3.2 Generalization ability

Figure 6 summarizes error statistics for the generalization tests described in
Section 2.3. (Detailed results for each test are shown in Figures 12, 13, 14,
15, and 16). The numerical errors in generalization testing datasets are
similar to the error in the training dataset in most cases, but outliers exist
in the generalization tests with errors substantially larger than that found
with the training dataset, and numerical instability does occur in some cases.
In most of the unstable cases, we could resolve the instability by using model
parameters from a training epoch different than the one that is optimal for
the training dataset; those cases are shown in bold letters in Figures 12, 13,
14, 15, and 16. The goal of the analysis here is to explore the performance
and flexibility of our method when trained on a small dataset, but we would
expect to see better generalization if we used a more diverse training
dataset.

Figure 6: Box plots of error statistics in generalization tests for (A) MAE,
(B) RMSE, and (C) R2. Each data point represents one combination of spatial
and temporal resolutions. Unstable simulations are not shown here but are
shown in Figures 12–16.

In Figure 6, the “3 times longer” case shows that predictive performance does
not significantly degrade over an extended prediction horizon. It also
suggests that the learned solver would work well in integrating wind fields
with similar characteristics (e.g., season and latitude) to the training
dataset.

The Dirac-delta initial condition test shows that changing the initial
condition deteriorates the performance when the new initial condition is
extraordinarily stiff (Figure 6). The normalized error is small, but this is
because the initial mass is small, rather than indicating good model
performance. The $R^{2}$ value better captures the large uncertainty of the
learned emulation in this case. In contrast, the learned solver is robust to
the initial condition with a smooth gradient (Gaussian shape; Figure 6).
Zhuang et al. (2021) found their learned 1-D advection scheme had a tendency
to have a square-shaped scalar even though they fed the Gaussian distribution-
like initial condition because their scheme was trained with the square wave
moving to one direction. In our study, that does not happen because our solver
could learn from not only the initial shape but also different spatial
patterns with time evolution. However, our results still imply that training
our model with a wider range of initial conditions would be useful.

The effects of seasonality and different spatial domains are larger than other
testing scenarios in terms of median error (Figure 6). The largest discrepancy
in median appears in changing the spatial domain to a longitudinal direction.
The largest outlier is seemingly shown in the October wind case, but this is
only because the Figure 6 does not show unstable simulations in the July and
$29.5\text{\,}\mathrm{\SIUnitSymbolDegree}$N tests for (2$\Delta$x,
64$\Delta$t) resolution. As above, we expect that these issues could be
resolved by training our model on a more diverse dataset that has different
ranges of velocity, stiffness of gradient, etc.

To look closer into the effect of wind velocity range on the generalization
skill of learned solvers, we plot MAE vs. maximum CFL number of the velocity
fields used in the generalization testing (Figure 7). This plot shows that the
learned solver tended to be unstable when the maximum CFL is far from the
maximum CFL of the training dataset. The testing cases with the largest
discrepancy in maximum CFL with the training dataset show numerical
instabilities in (2$\Delta$x, 64$\Delta$t) resolution. This feature might
originate from the structure of our learned solver. As shown in Figure 2, the
learned solver could scale the magnitude of scalar values but not in the case
of the velocity field. Therefore, the neural network could fail to give the
learned coefficients in a reasonable range when the input velocity is not
close to the training regime. Based on this intuition, in future work, we will
want to train the solver with a larger dataset that has a wider range of
velocity. However, as mentioned above, our goal here is to explore the
potential of this approach rather than producing a highly-tuned 1-D advection
solver.

Figure 7: Dependence of generalization test accuracy on the maximum CFL number
of target wind fields. Each box plot is the summary of errors in each testing
scenario integrating wind data different from the training dataset.

Overall, the results from our 1-D generalization tests suggest that our
learned scheme can work in a relatively wide set of circumstances, even when
trained on an extremely limited dataset. Our results additionally suggest that
it is important for the training dataset to include a distribution of CFL
numbers similar to the conditions the model is expected to be deployed under.
This will likely be even more crucial in 3-D applications where wind speed
varies greatly with altitude.

###  3.3 2-D demonstration

Figure 8 shows the performance of the learned 1-D advection solvers when
implemented in a 2-D simulation. Even though our learned solvers are trained
with 1-D, ten among the 35 2-D simulation cases achieve R2 values above 0.8.
This is a fascinating result given that the velocity fields used in 2-D
advection have more complex and diverse features than exist in our 1-D
training dataset. Simply, the 2-D advection requires 224 $\times$ 201 times of
1-D advection, and the relative variety of wind patterns in the 2-D dataset
could be a function of that dimension. In this regard, our approach looks
promising even in 2-D applications. The best accuracy is observed in (2$\Delta
x$, 8$\Delta t$), which shows remarkable performance including R2 = 1.00. The
maximal acceleration is 340 times in the coarsest resolution. Computational
time statistics for 2-D advection by the learned solvers are given in Table 5.

Figure 8: Performance of the learned solver in 2-D implementation in different
coarse-graining resolutions. Cases with bold numbers were unstable with the
optimal training epoch but stable when using a different training epoch.

Although the fact that our 1-D solver works at all in a 2-D context has
exceeded our expectations, we nonetheless proceed with an analysis of
strengths and weaknesses with the goal of informing future work. Unlike the
results in 1-D (Figure 3), 2-D simulations in several resolutions show poor
results (R${}^{2}<0.5$). There are unstable simulations in 2-D as well:
simulations in (1$\Delta x$, 2$\Delta t$), (1$\Delta x$, 4$\Delta t$),
(1$\Delta x$, 32$\Delta t$), and (2$\Delta x$, 64$\Delta t$) are unstable in
2-D advection, though they are stable in 1-D. We suspect that this limited
performance may be caused by the existence of conditions (for example a
greater range of wind speeds) that were not present in the training dataset.

To fairly evaluate our approach, we also test the 2-D simulations with the
low-resolution reference solver. Similar to 1-D coarsening cases, 2-D
advection using the reference solver shows good performance in the diagonal
elements of the coarsening performance matrix (Figure 9). In other coarsening
regimes, the reference solvers show a poor R2 but usually maintain small
errors. This is because the reference solver has the local slope limiter to
satisfy the total variance diminishing property, thus the simulation is free
from spurious oscillation when the simulation does not violate the CFL
condition (Durran, 2010). Instead, the solution dissipates to zero due to
numerical diffusion, resulting in low error but also low correlation.
Instability occurs when the coarsening results in the CFL number greater than
one. Numerical techniques exist that can handle conditions with CFL $>1$, but
they lie outside of the focus of the current study. The low-resolution
reference solvers also achieve very good performance in some resolutions
including (2$\Delta x$, 4$\Delta t$) and (2$\Delta x$, 8$\Delta t$). The
computational speedup by coarsening is always greater with the reference
solver than with the learned solver in the same resolution, as explained in
Section 3.1. The maximum acceleration in 2-D by the reference solver is
$7200\times$. Table 6 summarizes the computational time for 2-D advection by
the reference solvers in different scales.

Figure 9: Performance of the reference solver in 2-D at different resolutions.

We visualize the 2-D implementation results of both reference and learned
solvers at two resolutions (Figure 10). We select (2$\Delta x$, 8$\Delta t$)
and (16$\Delta x$, 64$\Delta t$), representing the best performance and the
maximum acceleration by the learned solvers respectively. In (2$\Delta x$,
8$\Delta t$), both solvers almost perfectly emulate the finest grid
simulation. Both solvers have a tendency to exaggerate the thickness of the
plume, but overall they represent the spatial pattern well. For (16$\Delta x$,
64$\Delta t$), we additionally feed the 1-D generalization testing datasets to
the model as well as the original training dataset to enhance the learned
model. Originally, the learned scheme trained only with the training dataset
have limited performance with R2 = 0.27 (Figure 8). With the increased
training data, the learned scheme R2 is doubled, showing the benefit of
introducing more data into the learned model. However, the learned scheme in
this resolution overestimates the concentration of the plume and sometimes
shows a negative concentration (with a magnitude less than $1$ ppb). These
artifacts result in a relatively large MAE value. The reference solver in this
resolution showed a similar range of R2, but dissipates the plume too quickly
and the overall concentration range near the plume is much too low. However,
even though the reference solver is dissipative at this resolution, the MAE
value was still smaller than the learned solver.

Figure 10: 2-D demonstration of coarse-graining by the learned and reference
advection schemes. The first row shows the results of the high-resolution
baseline solver. The second and fourth rows are the learned advection in
different resolutions. The third and fifth rows are the low-resolution
reference solvers at different resolutions.

##  4 Discussion

We will discuss the approach described herein in terms of strengths,
weaknesses, and directions for future research. The first noticeable strength
is the decent performance of the learned model in emulating the reference
model, even without using a large training dataset or hyperparameter tuning.
The coarse-grained learned solver achieves accuracy that is similar to or
better than the corresponding-resolution reference solver for conditions
similar to the training dataset and performs well in some cases where the
reference solver does not perform well, particularly for simulations with
temporal coarsening but without corresponding spatial coarsening. This is a
key benefit of our approach because for many applications of air quality
modeling, high spatial resolution is much more important than high temporal
resolution. For example, CTMs typically save results to disk every one or
three simulation hours, discarding all information at a higher temporal
resolution than that. Therefore, an approach that allows temporal coarse-
graining without corresponding spatial coarse-graining could achieve
computational acceleration without requiring any reduction in the detail of
the CTM output.

Furthermore, our learned solver generalizes fairly well to conditions it was
not trained on, even though it was trained on an extremely limited amount of
data. This good generalization performance extends to the application of our
learned model for 2-D simulations, even though it was only trained in 1-D.
Finally, our learned solvers can achieve substantial computational
acceleration relative to the reference solver, with maximum acceleration at
$\sim 10\times$ for 1-D and $\sim 100\times$ for 2-D simulations. Together,
these results suggest that a scaled-up version of the approach shown here
could be a good candidate for embedding within a full CTM for further testing,
and this is a priority area for future work.

One limitation of our learned advection model is that it is numerically
unstable under certain conditions. The number of scenarios where numerical
instability occurred was similar for the learned and reference solvers, but
the conditions that cause instability in the learned solver are somewhat
unpredictable, whereas with the reference solver instability predictably
occurs at conditions with the highest CFL values.

This limitation could be addressed in future work by using a larger and more
diverse training dataset to reduce the chance of the solver encountering
conditions that were not well represented during training, or by minimizing
error over a larger number of timesteps to encourage stability over longer
time horizons. Fine-tuning the hyperparameters could also help to develop a
more accurate and stable model. Other possibilities could be to bound in the
magnitude of the learned coefficients, similar to the idea of a flux limiter,
or to add terms to the training loss function to encourage mass conservation
or stability (e.g., the amplification ratio in von Neumann analysis should be
less than or equal to one).

We consider a learned solver to be useful when it is faster than the high-
resolution reference solver and more accurate than the equivalent-resolution
reference solver. A second limitation of the current study is that this is
only true for conditions where the reference solver is particularly
inaccurate, i.e. simulations with large timesteps and small grid cells or
simulations with small grid cells and large timesteps (Figure 4).

It may not ever be reasonable to expect a learned solver to be faster than the
reference solver when run at the same resolution because in general, a learned
model will have more degrees of freedom than a theoretically derived model,
and the additional degrees of freedom require additional computational
operations. However, future work could reduce the computational intensity of
the learned model by exploring different machine-learning architectures for
predicting the operator coefficients ($a$ and $b$ in Equation 8). A learned
model with a smaller per-step computational penalty as compared to the
reference solver would have a larger range of conditions where it was useful.

##  5 Conclusions

In this study, we train a learned 1-D advection solver with 10 days of
transport data and characterize its performance with the goal of exploring its
general strengths and weaknesses rather than optimizing its performance for a
particular use case. Even with this limited training dataset, the learned 1-D
advection scheme shows generally stable and accurate emulation which is robust
against unforeseen conditions. The 1-D learned solver could run the simulation
up to 18 times faster than the baseline solver and could be used in the 2-D
advection with an appropriate splitting technique. The 2-D implementation
results are still mostly stable, reasonably accurate, and overall more
successful than we expected ahead of time. Our approach demonstrates up to
340$\times$ computational speedup in the 2-D application. For simulations with
temporal coarse-graining but without spatial coarse-graining, numerical
instability occurred occasionally. Addressing this instability in future
research will help to fully realize the potential of this approach. We hope
that the findings in this study serve as an inspiration and a benchmark for
future studies on multidimensional learned advection operators.

###### Acknowledgements.

This work is supported by an Early Career Faculty grant from the National
Aeronautics and Space Administration (grant no. 80NSSC21K1813) and Assistance
Agreement RD-84001201-0 awarded by the U.S. Environmental Protection Agency.
It has not been formally reviewed by EPA. The views expressed in this document
are solely those of authors and do not necessarily reflect those of the
Agency. EPA does not endorse any products or commercial services mentioned in
this publication. MP is supported by the Carver Fellowship and Illinois
Distinguished Fellowship. The code to implement the reference and learned
advection and the velocity dataset are available at
https://doi.org/10.5281/zenodo.8111623. The output of the model training
including the model parameters and outputs will be available at
https://doi.org/10.13012/B2IDB-4743181_V1 and currently can be accessed via
https://databank.illinois.edu/datasets/IDB-4743181?code=8zCTeEEae1cO3_3JNLKndtCk4h9FwhzG2L6iCqhGzDU.
Appendix Tables and Figures

###  .1 Appendix Tables

Table 1: The best training epoch in each resolution tested in this study. In table cells with bold-face numbers, the training was not successful with 3-layer CNN but was successful with 4-layer CNN. Cases marked “unstable” were unable to perform stable simulations even with a 5-layer CNN. | $\Delta t$ | 2$\Delta t$ | 4$\Delta t$ | 8$\Delta t$ | 16$\Delta t$ | 32$\Delta t$ | 64$\Delta t$
---|---|---|---|---|---|---|---
$\Delta x$ | 95 | 112 | 13 | Unstable | Unstable | 225 | Unstable
2$\Delta x$ | 71 | 59 | 31 | 164 | 137 | 282 | 221
4$\Delta x$ | 57 | 33 | 75 | 197 | 100 | 199 | 100
8$\Delta x$ | 193 | 443 | 362 | 23 | 43 | 44 | 263
16$\Delta x$ | 65 | 166 | 85 | 31 | 41 | 25 | 56
|  |  |  |  |  |  |

Table 2: Computational time taken to integrate 10-day-long 1-D advection by
the learned advection solvers. (Unit: milliseconds; *: measured only once
because the computational time is large enough; N/A: not measured because the
solver is unstable)

| $\Delta t$ | 2$\Delta t$ | 4$\Delta t$ | 8$\Delta t$ | 16$\Delta t$ | 32$\Delta t$ | 64$\Delta t$
---|---|---|---|---|---|---|---
$\Delta x$ | 720.9 $\pm$ 9.8 | 365.5 $\pm$ 5.6 | 15722∗ | N/A | N/A | 1699 $\pm$ 83.3 | N/A
2$\Delta x$ | 407.3 $\pm$ 7.4 | 207.4 $\pm$ 8.6 | 104.6 $\pm$ 2.8 | 50.9 $\pm$ 2.3 | 27.7 $\pm$ 2.2 | 12.3 $\pm$ 1.4 | 6.6 $\pm$ 1.1
4$\Delta x$ | 252.2 $\pm$ 7.1 | 125.8 $\pm$ 2.8 | 62.5 $\pm$ 3.0 | 31.9 $\pm$ 2.9 | 15.9 $\pm$ 2.3 | 8.4 $\pm$ 1.6 | 3.8 $\pm$ 1.1
8$\Delta x$ | 173.0 $\pm$ 6.2 | 86.6 $\pm$ 4.2 | 43.1 $\pm$ 3.9 | 21.5 $\pm$ 2.9 | 10.7 $\pm$ 2.1 | 5.4 $\pm$ 1.5 | 2.8 $\pm$ 1.0
16$\Delta x$ | 134.6 $\pm$ 5.4 | 68.4 $\pm$ 7.7 | 33.7 $\pm$ 5.3 | 16.7 $\pm$ 3.8 | 8.3 $\pm$ 2.5 | 4.1 $\pm$ 1.7 | 2.1 $\pm$ 1.2
|  |  |  |  |  |  |

Table 3: Computational time taken to integrate 10-day-long 1-D advection by the learned advection solvers. (Unit: milliseconds; *: measured only once because the computational time is large enough; N/A: not measured because the solver is unstable) | $\Delta t$ | 2$\Delta t$ | 4$\Delta t$ | 8$\Delta t$ | 16$\Delta t$ | 32$\Delta t$ | 64$\Delta t$
---|---|---|---|---|---|---|---
$\Delta x$ | 37.9 $\pm$ 7.4 | 19.5 $\pm$ 6.1 | 9.5 $\pm$ 3.9 | 4.8 $\pm$ 2.8 | 2.4 $\pm$ 1.9 | 1.2 $\pm$ 1.4 | 0.6 $\pm$ 1.0
2$\Delta x$ | 19.9 $\pm$ 5.4 | 9.9 $\pm$ 3.9 | 5.0 $\pm$ 2.8 | 2.5 $\pm$ 2.0 | 1.2 $\pm$ 1.4 | 0.6 $\pm$ 1.0 | 0.3 $\pm$ 0.7
4$\Delta x$ | 11.2 $\pm$ 4.4 | 5.6 $\pm$ 3.1 | 2.8 $\pm$ 2.1 | 1.4 $\pm$ 1.6 | 0.7 $\pm$ 1.1 | 0.4 $\pm$ 0.8 | 0.2 $\pm$ 0.5
8$\Delta x$ | 7.0 $\pm$ 3.4 | 3.6 $\pm$ 2.6 | 1.8 $\pm$ 1.7 | 0.9 $\pm$ 1.2 | 0.4 $\pm$ 0.9 | 0.2 $\pm$ 0.6 | 0.1 $\pm$ 0.4
16$\Delta x$ | 4.8 $\pm$ 3.1 | 2.4 $\pm$ 2.1 | 1.2 $\pm$ 1.6 | 0.6 $\pm$ 1.0 | 0.3 $\pm$ 0.7 | 0.2 $\pm$ 0.5 | 0.1 $\pm$ 0.3
|  |  |  |  |  |  |

Table 4: The maximum CFL number of each wind field tested in this study. | $\Delta t$ | 2$\Delta t$ | 4$\Delta t$ | 8$\Delta t$ | 16$\Delta t$ | 32$\Delta t$ | 64$\Delta t$
---|---|---|---|---|---|---|---
$\Delta x$ | 0.23 | 0.45 | 0.91 | 1.82 | 3.63 | 7.26 | 14.52
2$\Delta x$ | 0.11 | 0.23 | 0.45 | 0.91 | 1.81 | 3.63 | 7.26
4$\Delta x$ | 0.06 | 0.11 | 0.23 | 0.45 | 0.90 | 1.81 | 3.61
8$\Delta x$ | 0.03 | 0.06 | 0.11 | 0.22 | 0.45 | 0.90 | 1.79
16$\Delta x$ | 0.01 | 0.03 | 0.05 | 0.11 | 0.22 | 0.43 | 0.87
|  |  |  |  |  |  |

Table 5: Computational time taken to integrate 10-day-long 2-D advection by
the learned advection solvers. (Unit: milliseconds; *: measured only once
because the computational time is large enough; N/A: not measured because the
solver is unstable)

| $\Delta t$ | 2$\Delta t$ | 4$\Delta t$ | 8$\Delta t$ | 16$\Delta t$ | 32$\Delta t$ | 64$\Delta t$
---|---|---|---|---|---|---|---
$\Delta x$ | 593,080∗ | 289,616∗ | 9,362,782∗ | N/A | N/A | N/A | N/A
2$\Delta x$ | 152,659∗ | 76,292∗ | 38,424∗ | 21,179∗ | 10,158∗ | 5,311∗ | N/A
4$\Delta x$ | 47,217∗ | 24,102∗ | 11,722∗ | 5,903∗ | 3,013 $\pm$ 59.7 | 1,576 $\pm$ 14.7 | 857.9 $\pm$ 6.8
8$\Delta x$ | 17,075∗ | 8,559∗ | 4,531 $\pm$ 405.4 | 2,124 $\pm$ 26.2 | 1,056 $\pm$ 12.9 | 530.5 $\pm$ 7.3 | 269.2 $\pm$ 7.8
16$\Delta x$ | 6,942∗ | 3,494 $\pm$ 31.3 | 1,740 $\pm$ 9.8 | 881.8 $\pm$ 11.2 | 439.6 $\pm$ 12.3 | 218.5 $\pm$ 8.5 | 110.4 $\pm$ 10.2
|  |  |  |  |  |  |

Table 6: Computational time taken to integrate 10-day-long 2-D advection by
the learned advection solvers. (Unit: milliseconds)

| $\Delta t$ | 2$\Delta t$ | 4$\Delta t$ | 8$\Delta t$ | 16$\Delta t$ | 32$\Delta t$ | 64$\Delta t$
---|---|---|---|---|---|---|---
$\Delta x$ | 37,442∗ | 19,531∗ | 9,162∗ | N/A | N/A | 1,323 $\pm$ 181.1 | N/A
2$\Delta x$ | 8,997∗ | 4,497 $\pm$ 42.2 | 2,391 $\pm$ 220.8 | 1,173 $\pm$ 9.1 | 640.6 $\pm$ 94.0 | 339.5 $\pm$ 5.8 | N/A
4$\Delta x$ | 2,485 $\pm$ 31.1 | 1,255 $\pm$ 7.5 | 647.8 $\pm$ 9.4 | 323.3 $\pm$ 9.8 | 167.8 $\pm$ 8.9 | 85.2 $\pm$ 2.1 | 44.2 $\pm$ 8.8
8$\Delta x$ | 826.1 $\pm$ 11.6 | 408.0 $\pm$ 7.1 | 203.1 $\pm$ 8.0 | 101.4 $\pm$ 9.0 | 50.8 $\pm$ 2.6 | 25.4 $\pm$ 6.0 | 12.6 $\pm$ 3.0
16$\Delta x$ | 329.8 $\pm$ 8.1 | 166.0 $\pm$ 5.8 | 83.4 $\pm$ 8.7 | 41.2 $\pm$ 7.7 | 21.0 $\pm$ 5.7 | 10.2 $\pm$ 3.6 | 5.2 $\pm$ 2.2
|  |  |  |  |  |  |

###  .2 Appendix Figures

Figure 11: Illustration of the mass-conserving downsampling method.

Figure 12: Summary of generalization testing with a longer time span. The
cases with bold numbers were unstable when using the optimal epoch for the
training data but stable when using a different epoch.

Figure 13: Summary of generalization testing with different initial
conditions. The cases with the bold numbers were unstable when using the
optimal epoch for the training data but stable when using a different training
epoch.

Figure 14: Summary of generalization testing with different latitudes. The
cases with the bold numbers were unstable when using the optimal epoch for the
training data but stable when using a different training epoch.

Figure 15: Summary of generalization testing with different seasons. The cases
with the bold numbers were unstable when using the optimal epoch for the
training data but stable when using a different training epoch.

Figure 16: Summary of generalization testing with the longitudinal
application. The cases with the bold numbers were unstable when using the
optimal epoch for the training data but stable when using a different training
epoch.

## References

  * Bey et al. (2001) Bey, I., and Coauthors, 2001: Global modeling of tropospheric chemistry with assimilated meteorology: Model description and evaluation. J. Geophys. Res. Atmos., 106 (D19), 23 073–23 095.
  * Bezanson et al. (2012) Bezanson, J., S. Karpinski, V. B. Shah, and A. Edelman, 2012: Julia: A fast dynamic language for technical computing. arXiv preprint arXiv:1209.5145.
  * Brunton et al. (2020) Brunton, S. L., B. R. Noack, and P. Koumoutsakos, 2020: Machine learning for fluid mechanics. Annu. Rev. Fluid Mech., 52, 477–508.
  * Courant et al. (1928) Courant, R., K. Friedrichs, and H. Lewy, 1928: Über die partiellen differenzengleichungen der mathematischen physik. Math. Ann., 100 (1), 32–74.
  * Durran (2010) Durran, D. R., 2010: Numerical methods for fluid dynamics: With applications to geophysics, Vol. 32. Springer Science & Business Media.
  * Eastham et al. (2018) Eastham, S. D., and Coauthors, 2018: Geos-chem high performance (gchp v11-02c): a next-generation implementation of the geos-chem chemical transport model for massively parallel applications. Geosci. Model Dev., 11 (7), 2941–2953.
  * EPA (2022) EPA, 2022: Cmaq. Zenodo, 10.5281/zenodo.7218076.
  * Grell et al. (2005) Grell, G. A., S. E. Peckham, R. Schmitz, S. A. McKeen, G. Frost, W. C. Skamarock, and B. Eder, 2005: Fully coupled “online” chemistry within the wrf model. Atmos. Environ., 39 (37), 6957–6975.
  * Hendrycks and Gimpel (2016) Hendrycks, D., and K. Gimpel, 2016: Gaussian error linear units (gelus). arXiv preprint arXiv:1606.08415.
  * Huang and Seinfeld (2022) Huang, Y., and J. H. Seinfeld, 2022: A neural network-assisted euler integrator for stiff kinetics in atmospheric chemistry. Environ. Sci. Technol., 56 (7), 4676–4685.
  * Innes et al. (2019) Innes, M., A. Edelman, K. Fischer, C. Rackauckas, E. Saba, V. B. Shah, and W. Tebbutt, 2019: A differentiable programming system to bridge machine learning and scientific computing. arXiv preprint arXiv:1907.07587.
  * Innes et al. (2018) Innes, M., and Coauthors, 2018: Fashionable modelling with flux. arXiv preprint arXiv:1811.01457.
  * JuliaCI (2022) JuliaCI, 2022: BenchmarkTools.jl. URL https://github.com/JuliaCI/BenchmarkTools.jl, accessed: 2023.07.27.
  * Keller and Evans (2019) Keller, C. A., and M. J. Evans, 2019: Application of random forest regression to the calculation of gas-phase chemistry within the geos-chem chemistry model v10. Geosci. Model Dev., 12 (3), 1209–1225.
  * Keller et al. (2017) Keller, C. A., M. J. Evans, J. N. Kutz, and S. Pawson, 2017: Machine learning and air quality modeling. 2017 IEEE International Conference on Big Data (Big Data), IEEE, 4570–4576.
  * Kelp et al. (2020) Kelp, M. M., D. J. Jacob, J. N. Kutz, J. D. Marshall, and C. W. Tessum, 2020: Toward stable, general machine-learned models of the atmospheric chemical system. J. Geophys. Res. Atmos., 125 (23), e2020JD032 759.
  * Kelp et al. (2022) Kelp, M. M., D. J. Jacob, H. Lin, and M. P. Sulprizio, 2022: An online-learned neural network chemical solver for stable long-term global simulations of atmospheric chemistry. J. Adv. Model. Earth Syst., 14 (6), e2021MS002 926.
  * Kim et al. (2021) Kim, S., W. Ji, S. Deng, Y. Ma, and C. Rackauckas, 2021: Stiff neural ordinary differential equations. Chaos, 31 (9), 093 122.
  * Kingma and Ba (2014) Kingma, D. P., and J. Ba, 2014: Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980.
  * Kochkov et al. (2021) Kochkov, D., J. A. Smith, A. Alieva, Q. Wang, M. P. Brenner, and S. Hoyer, 2021: Machine learning–accelerated computational fluid dynamics. Proc. Natl. Acad. Sci. (USA), 118 (21), e2101784 118.
  * Lin et al. (1994) Lin, S.-J., W. C. Chao, Y. Sud, and G. Walker, 1994: A class of the van leer-type transport schemes and its application to the moisture transport in a general circulation model. Mon. Wea. Rev., 122 (7), 1575–1593.
  * Lin and Rood (1996) Lin, S.-J., and R. B. Rood, 1996: Multidimensional flux-form semi-lagrangian transport schemes. Mon. Wea. Rev., 124 (9), 2046–2070.
  * NASA (2022) NASA, 2022: GMAO Product GEOS Near-Real Time Data Products. URL https://gmao.gsfc.nasa.gov/GMAO˙products/NRT˙products.php, accessed: 2023.07.27.
  * Philip et al. (2016) Philip, S., R. V. Martin, and C. A. Keller, 2016: Sensitivity of chemistry-transport model simulations to the duration of chemical and transport operators: a case study with geos-chem v10-01. Geosci. Model Dev., 9 (5), 1683–1695.
  * Shen et al. (2022) Shen, L., D. J. Jacob, M. Santillana, K. Bates, J. Zhuang, and W. Chen, 2022: A machine-learning-guided adaptive algorithm to reduce the computational cost of integrating kinetics in global atmospheric chemistry models: application to geos-chem versions 12.0. 0 and 12.9. 1. Geosci. Model Dev., 15 (4), 1677–1687.
  * Stachenfeld et al. (2021) Stachenfeld, K., and Coauthors, 2021: Learned coarse models for efficient turbulence simulation. arXiv preprint arXiv:2112.15275.
  * Sturm et al. (2023) Sturm, P. O., A. Manders, R. Janssen, A. Segers, A. S. Wexler, and H. X. Lin, 2023: Advecting superspecies: Efficiently modeling transport of organic aerosol with a mass-conserving dimensionality reduction method. J. Adv. Model. Earth Syst., 15 (3), e2022MS003 235.
  * Vinuesa and Brunton (2022) Vinuesa, R., and S. L. Brunton, 2022: Enhancing computational fluid dynamics with machine learning. Nat. Comput. Sci., 2 (6), 358–366.
  * Zhuang et al. (2020) Zhuang, J., D. J. Jacob, H. Lin, E. W. Lundgren, R. M. Yantosca, J. F. Gaya, M. P. Sulprizio, and S. D. Eastham, 2020: Enabling high-performance cloud computing for earth science modeling on over a thousand cores: application to the geos-chem atmospheric chemistry model. J. Adv. Model. Earth Syst., 12 (5), e2020MS002 064.
  * Zhuang et al. (2021) Zhuang, J., D. Kochkov, Y. Bar-Sinai, M. P. Brenner, and S. Hoyer, 2021: Learned discretizations for passive scalar advection in a two-dimensional turbulent flow. Phys. Rev. Fluids, 6 (6), 064 605.


# Individualized conformal

Fernando Delbianco - Fernando Tohmé
Departamento de Economía  Universidad Nacional del Sur
Instituto de Matemática de Bahía Blanca  UNS - CONICET

###### Abstract

The problem of individualized prediction can be addressed using variants of
conformal prediction, obtaining the intervals to which the actual values of
the variables of interest belong. Here we present a method based on detecting
the observations that may be relevant for a given question and then using
simulated controls to yield the intervals for the predicted values. This
method is shown to be adaptive and able to detect the presence of latent
relevant variables.
Keywords: Conformal Prediction, Individualized Inference, Split and Jacknife
Distribution-Free Inference.

##  1 Introduction

We present a novel approach to individualized conformal prediction, intended
to yield the intervals in which the actual predicted values of certain
variables might be found. We proceed by first detecting relevant observations
for a given target query using a technique reminiscent to divide and conquer
(Chen et al., (2021)). Then we apply a data augmentation procedure similar to
the generation of repro samples (Xie and Wang, (2022)), simulating controls as
if they were bootstrapped (as in Tran et al., (2017)). The development of this
approach is based on the statistical technique known as Conformal Prediction.

Conformal Prediction (CP) is a distribution-free prediction methodology
usually based on machine learning systems. CP generates predictions about new
test data points on the basis of training labeled datasets, only assuming the
exchangeability of the data. It requires a pre-specified level that restricts
the frequency of errors that the algorithm is allowed to make (see, for
instance, Vovk et al., (2022), Xie and Zheng, (2022), Lei et al., (2018),
Barber et al., (2021)). Variants of this methodology are split and jacknife
distribution-free conformal inference (Lei et al., (2018)). A recent
development that is closely related to our contribution is the adaptive
version of CP, based on the application of self-supervised learning (Seedat et
al., (2023)).

Another strand in the literature that is relevant for our approach is
individualized inference. This is a particularly useful methodology in the
case of large datasets, since it addresses the heterogeneity of data. One
approach is based on generating iGroups of data sharing common features with
targeted individuals (Cai et al., (2021)). Alternatively, in the iFusion
approach the group of relevant data is generated by fusing the inferences from
individuals that are similar to the targeted one (Shen et al., (2020)).
Individualized inference is also invoked to categorize cases by mixing
Gaussian processes generated by individuals with related features (Alaa and
Van Der Schaar, (2017)).

The specific techniques used to carry out the inferences are varied.
Chernozhukov et al., (2017) develops a Debiased Machine Learning method to
obtain inferences about specific parameters. More traditional methods can also
be applied, like inferring confidence (Xie and Singh, (2013), Schweder and
Hjort, (2016)) and predictive (Shen et al., (2018)) distributions, which could
be seen as “Bayesian-like frequentist techniques” that derive distributions up
from observed data. In practical terms, an important contribution is the
implementation in the R language of trainable $p$-value functions (instead of
just $p$-values) by Infanger and Schmidt-Trucksäss, (2019).

Our own approach has the following features:

  * •

It is adaptive.

  * •

It admits new unlabeled data as input.

  * •

It satisfies the condition of exchangeability by restricting the focus on
relevant, and then simulated, data.

  * •

It detects empirically the presence of latent relevant variables, by
implicitly getting rid of the factor that generates a failure of the i.i.d.
condition.

This paper is structured as follows. We first present, in Section 2, the
motivation of our proposal. Section 3 details a procedure implementing the
method developed. Section 4 presents the empirical setting in which we test
the methodology and the corresponding results. Section 5 concludes by
presenting ideas for further developments.

##  2 Motivation

We assume, as in Delbianco et al., (2021), a statistical model of a data
generation process, which can be described as $\\{\mathcal{O},\mathcal{P}\\}$,
where $\mathcal{O}$ is the set of observations while
$\mathcal{P}=\\{P_{\theta}:\theta\in\Theta\\}$ is a family of probability
distributions over $\mathcal{O}$ and $\Theta$ is the space of parameters of
the model. The goal is again to estimate intervals for the parameters in
response to any query q, where q is a specific request for information under a
given inference method $\mathcal{I}$ applied on the database.

The query defines several dimensions. First, $\mathcal{O}$ consists of entries
$\\{o_{i}=\langle x_{i},y_{i}\rangle\\}_{i=1,\ldots,n}$, where $x_{i}$ is a
matrix of $p$ variables and $n$ observations, the tail of $o_{i}$, while
$y_{i}$ is the head, a vector of dimension $n$. The query $q$ consists of a
tail, $x_{0}$, with no head, and with $x_{o}\neq x_{i}$, for $i=1,\ldots,n$.
We assume that there exists a class of latent variables $\mathcal{S}$ such
that given a query $q$ there exists a corresponding $s_{q}\in\mathcal{S}$
yields a class of relevant observations $\mathcal{O}_{q}\subseteq\mathcal{O}$
111Another reason to use controls from a subset of the observations is due the
notion of the Law of Large Populations. Meng, (2018) refers to the difference
between data quantity and data quality..

Figure 1: Latent relevant area $S_{q}$$y$$x$$q$ A particular value of a tail
$q$, is associated to a latent variable $s_{q}$ that yields the relevant space
of tails $x\in X_{q}$ and as a consequence allows to infer the corresponding
class of heads $y\in Y_{q}$. Then, this bi-dimensional figure can be extended
to a tri-dimensional one, with the sequence of different queries as a third
axis. This represents how the relevant set varies according to $q$, learning
the associated $\theta^{q}$.

We then generate a class of controls $\bar{\mathcal{O}}_{q}$ verifying that
$\mathcal{O}_{q}\subseteq\bar{\mathcal{O}}_{q}$. Based on
$\bar{\mathcal{O}}_{q}$ the application of the inference procedure
$\mathcal{I}$ yields an interval $\mathbf{I}_{q}$, with
$\theta^{q}\in\mathbf{I}_{q}$ such that $\bar{\mathcal{O}}_{q}$ can be
understood as a set of draws from a distribution
$\bar{P}_{\theta^{q}}\in\mathcal{P}$ 222Also known as Transitional inference,
as in Li and Meng, (2021).. The relevant set of observations is given by
$\mathcal{O}_{q}=f_{\mathcal{S}}(s_{q})$, where
$f_{\mathcal{S}}:\mathcal{S}\rightarrow 2^{\mathcal{O}}$ characterizes a
selection procedure.

Figure 2: Simulated controls
$\in\bar{\mathcal{O}}_{q}$$\mathcal{O}$$\mathcal{O}_{q}$$\bar{\mathcal{O}}_{q}$$\bar{P}_{\theta^{q}}$$P_{\theta^{q}}$
The relevant set $\mathcal{O}_{q}$, detected by means of $f_{\mathcal{S}}$, is
extended to an enlarged set, $\bar{\mathcal{O}}_{q}$. The difference in the
inference based in new controls can be measured as
$\Delta_{\theta^{q}}=E(\bar{P}_{\theta^{q}})-E(P_{\theta^{q}})$.

For an example consider the case of different economies $i$, each one
described by a vector of a few macroeconomic variables $o_{i}=\langle
x_{i},y_{i}\rangle$, where $y_{i}$ is the GDP of country $i$. We can ask, for
any given economy $\bar{o}\notin\mathcal{O}$ what is its expected GDP in five
years, knowing only $\bar{x}$. Now assume that a latent variable is the
productivity of the leading sectors in an economy. Then we can group all the
countries in $\mathcal{O}$ with a similar productivity as that of $\bar{o}$
and generate a class of controls based on this choice.

Taking another example, we can consider a dataset of programs, students and
grades. For a new student, we can ask different questions. For instance, how
long will take for a new student to finish her studies? What is the
probability of her changing majors? or the probability of dropping out
school?, among many other possible queries about this particular student that
could be answered with the observations in the dataset. Of course, a different
query may require different controls. And this can be true for both the
variables and the observations. So the $s_{q}$ latent variable will yield the
portion of observations that is relevant given $q$, and will let us estimate
$f_{\mathcal{S}}$ in order to simulate new controls and gain robustness in the
inference.

To proceed in this way we need to specify, for each $q$ a latent variable (or
set of latent variables) that are relevant for the query. Then, we have to
distinguish its range using some proxy or measurement on the available
$\mathcal{O}$. This will yield the class of relevant observations. For the
sake of simplicity, we could assume that exists only two types of query,
$q_{A}$ and $q_{B}$, and the latent variable $s_{q}$ defines
$\mathcal{O}_{q_{A}}$ and $\mathcal{O}_{q_{B}}$. This setting can be later
generalized to cases of more than two types of query. In the particular
example of the students, let us assume that there are two types of students,
associated to a latent variable (which can capture, for instance,
socioeconomic or cognitive advantages). Then, each $q$ will be associated to
its corresponding type of student and according to $f_{\mathcal{S}}$, mapped
to the class of observations of students of that type.

As said, the class $\mathcal{O}_{q}$ is that of the entries in $\mathcal{O}$
that are relevant to answering $q$. But once obtained these relevant
observations, a robust inference requires generating new controls, not present
in $\mathcal{O}_{q}$. This is achieved by creating pairs $\langle x,y\rangle$
similar to those in $\mathcal{O}_{q}$ but without assuming that they share
with the latter a common value of the latent variable. This means that only
the observable features of the entries $o\in\mathcal{O}_{q}$ must be used to
create fictitious controls in $\bar{\mathcal{O}}_{q}$, for instance as in the
examples of Section 5 of Delbianco et al., (2021).

##  3 Framework

As said, our procedure takes elements from different approaches:

  * •

Lei et al., (2018) presents three versions of the conformal prediction method.
One is based on the original formulation of Vovk et al., (2022). The other two
are the split conformal and the jackknife alternatives.

  * •

Lei et al., (2018) also presents an R code for distribution-free prediction,
which we modify to be applied in our project.

  * •

We use both a classical linear procedure and one based on elastic nets.
Alternatively, we can include a Gaussian Kernel regression to implement a non-
parametric approach. Other methods can be also easily applied using the
estimations and predictions generated by the R package Caret (Kuhn, (2008)),
which is convenient for the following additional reasons:

    1. 1.

It implements a previous training stage, cross-validating the results in order
to choose the meaningful features. It even allows to apply methodologies that
split the entries of the training database in a different way as conformal
prediction. This allows to make a finer selection of the $n_{r}$ relevant
entries and the $2n_{r}$ controls.

    2. 2.

It allows to choose from almost fifty prediction methodologies, among which
are GLM, Kernel, Quantile regression, Random Forest, etc. But, as said, we use
in a first run only regression and LASSO.

  * •

The parameters in our exercise are: $\alpha$ for confidence intervals, $\rho$
for splitting in split conformal prediction, and $\gamma$ for cosine
similarity (for percentile similarity we also use $\alpha$).

  * •

We add two stages of prediction, using the concept of relevance for individual
inference of our original presentation (Delbianco et al., (2021)):

    1. 1.

The first stage evaluates the relevance of entries in the database for each
new observation, based on the degree of similarity with the rest of the
database. The relevance is defined, when $p$ is small, in terms of the
percentage of distance between the tails of observations and those in the
database. Otherwise, when the dimension of the tails is large (in particular
when $p>n$), we apply a cosine-based measure. That is, the closer to 1 the
cosine between the respective tails is, the more relevant is an entry for a
new observation. A prediction can be made on the basis of the relevant
entries: a possible value of the head of the new observation can be predicted,
using a standard prediction method on the subset of relevant entries of the
original database.

    2. 2.

A second type of prediction is based on the generation of controls using the
relevant entries collected in the previous stage. For each individual entry in
a query for $X_{0}$, we generate another tail-head vector by adding a small
noise to its features. Thus, given $n_{r}$ relevant entries ($n_{r}<<n$), we
obtain $2n_{r}$ control vectors at this stage.

We will now present the pseudocode of the algorithm that yields the intervals
of prediction for the queries. Notice that for simplicity we do not state
explicitly the steps and parameters of the methods drawn from $\mathcal{M}$
(Lei et al., (2018), implemented in the R package
ConformalInference333https://github.com/ryantibs/conformal), $\mathcal{A}$
(implemented in the Caret package444https://topepo.github.io/caret/index.html)
and the similarity functions in $\mathcal{S}$.

Algorithm 1 Individualized Conformal Prediction

Input: $X$, $Y$, $X_{0}$, Regression method $\mathcal{A}$, Similarity method
$\mathcal{S}$, Conformal method $\mathcal{M}$, miscoverage level $\alpha$,
Similarity intensity $\gamma$

Output: Prediction intervals for each element in $X_{0}$

$\alpha\leftarrow\alpha_{i}\in\\{0;1\\}$

$\gamma\leftarrow\gamma_{i}\in\\{0;1\\}$

$\mathcal{A}\leftarrow\mathcal{A}_{i}$ $\triangleright$ OLS, LASSO, Kernel

$\mathcal{M}\leftarrow\mathcal{M}_{i}$ $\triangleright$ Conformal, Split,
Jackkinfe

$\mathcal{S}\leftarrow\mathcal{S}_{i}$ $\triangleright$ Percentile, Cosine

for each $x\in X_{0}$ do

1.1 Predict $y_{0}$ with $\mathcal{A}$ and $X$

1.2 Use residuals and $\mathcal{M}$ $\rightarrow$ interval ($\mathcal{C}$)

2.1 Use $\gamma$ and $\mathcal{S}$ $\rightarrow$ $X_{rel}$

2.2 Predict $y_{0}$ with $\mathcal{A}$ and $X_{rel}$

2.3 Use residuals and $\mathcal{M}$ $\rightarrow$ relevant interval
($\mathcal{C}_{rel}$)

3.1 Simulate $X_{sim}$ mimicking the distribution of $X_{rel}$ and
$n{sim}=n{rel}$

3.2 Predict $y_{0}$ with $\mathcal{A}$ and $X_{sim}$

3.3 Use residuals and $\mathcal{M}$ $\rightarrow$ relevant and simulated
interval ($\mathcal{C}_{sim}$)

end for each $x\in X_{0}$

Return: $\mathcal{C}$, $\mathcal{C}_{rel}$, $\mathcal{C}_{sim}$

A diagrammatic representation of this algorithm is depicted in Figure 3. We
can see that there are three paths, leading either to $\mathcal{C}$,
$\mathcal{C}_{rel}$ or $\mathcal{C}_{sim}$. Each of those can be executed
independently, although the latter two share a stage of detection of relevant
observations.

In Figure 3 the main decisions to be made are represented by the red circle.
Each single change in the choices made can lead to different results. In
principle we do not know how robust are the results to small differences in
the choices made, but some are easily predicted. So, for instance, a lower
significance level $\alpha$ yields longer intervals while a higher similarity
level $\gamma$ reduces the number of relevant observations. Other than that,
the actual choices of estimation method, type of conformal inference and
relevance criterion depend on the problems at hand and how the user assesses
them.

$D=(X,Y)$
---
$q=x_{0}$
$\mathcal{A}$
---
$\mathcal{S}$
$\mathcal{M}$
$\alpha$
$\gamma$
$\hat{y_{0}}=\mathcal{A}(X)$$\hat{\mathcal{C}}$$D_{rel}\in
D$$\hat{y_{0}}=\mathcal{A}(X_{rel})$$\hat{\mathcal{C}}_{rel}$$X_{sim}\sim
X_{rel}$$\hat{y_{0}}=\mathcal{A}(X_{sim})$$\hat{\mathcal{C}}_{sim}$ Figure 3:
Diagram of Algorithm 1. The red circle encloses the initialization stage of
the algorithm in which the estimation method, the type of conformal inference
and the relevance criterion are chosen. These decisions impact on all the
outputs. Up from this point, three alternative paths can be taken. The first
one does not involve relevant data and yields the standard result of conformal
prediction, $\mathcal{C}$. The second and the third paths require selecting a
subset of observations of $D$ that are similar to the query $x_{0}$. From then
on, two possible alternatives consist in either using just that subset to
obtain the confidence interval $\mathcal{C}_{rel}$ or to generate a synthetic
sample enlarging the class of relevant data to obtain $\mathcal{C}_{sim}$. In
Algorithm 1 these alternative paths are denoted 1, 2 and 3, with their
respective steps (1.1, 1.2, 2.1, etc.)

To compare the results of the different methods chosen we apply four metrics.
Two of them, $A$ and $B$, do not evaluate the intervals but the predicted
values of the variable of interest.

  * $A$.

Distance from the forecast: the absolute value of the difference between the
true value $y_{0}$ and the value predicted by $\mathcal{A}$ on $X$, $X_{rel}$,
or $X_{sim}$. Naturally, the smaller $A$ the more accurate the prediction.

  * $B$.

Distance from the forecast, as a percentage of $y_{0}$: it is defined as
$\frac{A}{y_{0}}$ indicating how far is the forecast from the actual one in
proportion to the latter. In this case it is not clear that smaller is always
better.

  * $C$.

Length of the interval: the difference between the upper and lower limits of
the intervals generated by Algorithm 1. The length of an interval depends on
$\alpha$, but at the same significance level, a shorter interval indicates a
lower uncertainty about the forecast.

  * $D$.

Normalized distance from the forecast: it is defined as $\frac{A}{C}$. It
indicates how large is the error with respect the length of the interval.

Two additional measures may also provide information about the quality of the
intervals generated by Algorithm 1, namely coverage and excess (see Seedat et
al., (2023)). In the exercises reported in Section 4, all the methods make
similar predictions either inside or outside the confidence interval, and thus
these metrics may not provide extra information. But for larger heterogeneous
databases they may become more informative.

##  4 Empirical setting

For our preliminary explorations we generate a dataset of 750 observations,
where each third corresponds to a different and heterogeneous data generating
process. The corresponding settings are presented in Table 1.

We also enlarge the number of variables in the tails of the entries in the
databases ($p$). The description of the models is shown in Table 2.

Tables 3 and 4 show the results of running Algorithm 1 on the simulated data
of Table 1. The metrics on these results as well as on those obtained applying
Gaussian Kernel and LASSO are presented in Tables 5, 6 and 7.

Figures 4, 5 and 6 show that the predicted intervals are adapted to the true
values. Notice that the residuals increase at the tails of the distribution of
true values. This is not quite surprising since we do not use $y_{0}$ to infer
the intervals of the new data used to make the final forecast.

(a) Predictions

(b) Residuals

Figure 4: Conformal prediction with simulated relevant controls and OLS

(a) Predictions

(b) Residuals

Figure 5: Conformal prediction with simulated relevant controls and LASSO

(a) Predictions

(b) Residuals

Figure 6: Conformal prediction with simulated relevant controls and Kernel
regression

Figure 7 is obtained using the smaller dataset to show the comparison among
the three conformal prediction methods. We can see that they do not differ
much for the three values used as queries.

(a) Predictions

(b) Residuals

Figure 7: Results of the three methods with simulated relevant controls and
LASSO

##  5 Discussion

As already noted, the intervals adapt to each query. The point forecasts
depend more on the data used than on the method of conformal prediction
chosen. The prediction of the point values tends to be very conservative and
thus the residuals become negative on the left and positive on the right. In
most cases the predicted intervals capture the true values.

The cases analyzed empirically here are, of course, very simple. They do not
allow us to explore exhaustively all the consequences of the choices at the
initialization phase of Algorithm 1. But we can make the preliminary
observation that no clear winner can be ascertained among the conformal
prediction method, since the three variants yield in average similar results,
although the split conformal one is less computationally costly.

Our method yields adaptive intervals using only relevant plus simulated data
as to ensure the robustness of the predictions. This divide and conquer
strategy is particularly appropriate for large and heterogeneous datasets.
Subsets of data individualize the forecasts of variables of which specific
information is not available.

As an extra bonus, for each query we obtain a distribution of relevant data.
This captures a latent variable (this facilitates the prediction of the
intervals corresponding to the query), providing useful information about the
individuals in the same “class”.

On the downside, we have seen that residuals are more widely dispersed at the
tails of the distribution of actual values. This is inherited from the
application of conformal methods, unlike in the self-supervised method of
Seedat et al., (2023), where residuals are used to train the system to improve
its predictions.

Further work is needed to make a more precise assessment of the pros and cons
of the method, enriching it by relating it to other procedures presented in
the literature. Possible lines of study are:

  * •

Run more experiments using a larger $p$ or a higher-dimensional setup.

  * •

Marginalize the individual result up from that of the the relevant group (or
even up from the group of non-relevant ones) as in Zhou et al., (2023).

  * •

Explore whether our individualization strategy is compatible with the version
of divide and conquer discussed in Chen et al., (2021).

  * •

Use residuals to create corrections at the tails of the intervals as to reduce
the epistemic uncertainty at the boundaries Alaa et al., (2023).

## References

  * Alaa et al., (2023) Alaa, A. M., Hussain, Z., and Sontag, D. (2023).  Conformalized unconditional quantile regression.  arXiv preprint arXiv:2304.01426.
  * Alaa and Van Der Schaar, (2017) Alaa, A. M. and Van Der Schaar, M. (2017).  Bayesian inference of individualized treatment effects using multi-task gaussian processes.  Advances in neural information processing systems, 30.
  * Barber et al., (2021) Barber, R. F., Candes, E. J., Ramdas, A., and Tibshirani, R. J. (2021).  Predictive inference with the jackknife+.  The Annals of Statistics, 49(1):486–507.
  * Cai et al., (2021) Cai, C., Chen, R., and Xie, M.-g. (2021).  Individualized group learning.  Journal of the American Statistical Association, pages 1–17.
  * Chen et al., (2021) Chen, X., Cheng, J. Q., and Xie, M.-g. (2021).  Divide-and-conquer methods for big data analysis.  arXiv preprint arXiv:2102.10771.
  * Chernozhukov et al., (2017) Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., and Newey, W. (2017).  Double/debiased/neyman machine learning of treatment effects.  American Economic Review, 107(5):261–65.
  * Delbianco et al., (2021) Delbianco, F., Fioriti, A., and Tohmé, F. (2021).  A methodology to answer to individual queries: finding relevant and robust controls.  Behaviormetrika, 48(2):259–282.
  * Infanger and Schmidt-Trucksäss, (2019) Infanger, D. and Schmidt-Trucksäss, A. (2019).  P value functions: An underused method to present research results and to promote quantitative reasoning.  Statistics in Medicine, 38(21):4189–4197.
  * Kuhn, (2008) Kuhn, M. (2008).  Building predictive models in r using the caret package.  Journal of statistical software, 28:1–26.
  * Lei et al., (2018) Lei, J., G’Sell, M., Rinaldo, A., Tibshirani, R. J., and Wasserman, L. (2018).  Distribution-free predictive inference for regression.  Journal of the American Statistical Association, 113(523):1094–1111.
  * Li and Meng, (2021) Li, X. and Meng, X.-L. (2021).  A multi-resolution theory for approximating infinite-p-zero-n: Transitional inference, individualized predictions, and a world without bias-variance tradeoff.  Journal of the American Statistical Association, 116(533):353–367.
  * Meng, (2018) Meng, X.-L. (2018).  Statistical paradises and paradoxes in big data (i) law of large populations, big data paradox, and the 2016 us presidential election.  The Annals of Applied Statistics, 12(2):685–726.
  * Schweder and Hjort, (2016) Schweder, T. and Hjort, N. L. (2016).  Confidence, likelihood, probability, volume 41.  Cambridge University Press.
  * Seedat et al., (2023) Seedat, N., Jeffares, A., Imrie, F., and van der Schaar, M. (2023).  Improving adaptive conformal prediction using self-supervised learning.  arXiv preprint arXiv:2302.12238.
  * Shen et al., (2018) Shen, J., Liu, R. Y., and Xie, M.-g. (2018).  Prediction with confidence—a general framework for predictive inference.  Journal of Statistical Planning and Inference, 195:126–140.
  * Shen et al., (2020) Shen, J., Liu, R. Y., and Xie, M.-g. (2020).  i fusion: Individualized fusion learning.  Journal of the American Statistical Association, 115(531):1251–1267.
  * Tran et al., (2017) Tran, T., Pham, T., Carneiro, G., Palmer, L., and Reid, I. (2017).  A bayesian data augmentation approach for learning deep models.  Advances in neural information processing systems, 30.
  * Vovk et al., (2022) Vovk, V., Gammerman, A., and Shafer, G. (2022).  Algorithmic learning in a random world.  Springer.
  * Xie and Singh, (2013) Xie, M.-g. and Singh, K. (2013).  Confidence distribution, the frequentist distribution estimator of a parameter: A review.  International Statistical Review, 81(1):3–39.
  * Xie and Wang, (2022) Xie, M.-g. and Wang, P. (2022).  Repro samples method for finite-and large-sample inferences.  arXiv preprint arXiv:2206.06421.
  * Xie and Zheng, (2022) Xie, M.-g. and Zheng, Z. (2022).  Homeostasis phenomenon in conformal prediction and predictive distribution functions.  International Journal of Approximate Reasoning, 141:131–145.
  * Zhou et al., (2023) Zhou, Z., Li, D., Huh, D., Xie, M., and Mun, E.-Y. (2023).  A simulation study of the performance of statistical models for count outcomes with excessive zeros: Focusing on the marginalized model.  arXiv preprint arXiv:2301.12674.

## Appendix

Setting $A$ |  Setting $B$ |  Setting $C$
---|---|---
one variable | two variables | heteroskedastic
$n_{1}=250$ | $n_{2}=250$ | $n_{3}=250$
$u_{1}=\mbox{rnorm}(n_{1})$ | $u_{2}=\mbox{rnorm}(n_{2})$ | $x_{1c}=\mbox{rnorm}(n_{3},1,1)$
$x_{1a}=\mbox{rnorm}(n_{1},1,1)$ | $x_{1b}=\mbox{rnorm}(n_{2},3,1)$ | $u_{3}=\mbox{rnorm}(n_{3})*\frac{x_{1c}}{2}$
$x_{2a}=\mbox{rnorm}(n_{1},2,1)$ (irrelevant) | $x_{2b}=\mbox{rnorm}(n_{2},2,2)$ | $x_{2c}=\mbox{rnorm}(n_{3},3,2)$
$y_{a}=0.5*x_{1a}+u_{1}$ | $y_{b}=0.5*x_{1b}+0.33*x_{2b}+u_{2}$ | $y_{c}=0.5*x_{1c}*x_{1c}+0.33*x_{2c}+u_{3}$
new data for $A$ |  new data for $B$ |  new data for $C$
$x_{01a}=\mbox{rnorm}(1,1,1)$ | $x_{01b}=\mbox{rnorm}(1,3,1)$ | $x_{01c}=\mbox{rnorm}(1,1,1)$
$x_{02a}=\mbox{rnorm}(1,2,2)$ | $x_{02b}=\mbox{rnorm}(1,2,2)$ | $x_{02c}=\mbox{rnorm}(1,3,2)$
$y_{0a}=0.5*x_{01a}$ | $y_{0b}=0.5*x_{01b}+0.33*x_{02b}$ | $y_{0c}=0.5*x_{01c}+0.33*x_{02c}$
$+$ $\mbox{rnorm}(1)$ |  $+$ $\mbox{rnorm}(1)$ |  $+$ $\mbox{rnorm}(1)*\frac{x_{01c}}{2}$
Table 1: Simulation datasets with large $n$ and small $p$ (using functions from the R language) $DGP_{1}$: |  $DGP_{2}$: |  $DGP_{3}$:
---|---|---
12 variables $N(0,1)$ with 2 coefficients $=0$ | shift in the mean of the features | shift in the mean of the coefficients
$n_{1}=100;p_{1}=12;s_{1}=2$ | $n_{2}=100;p_{2}=12;s_{2}=2$ | $n_{3}=100;p_{3}=12;s_{3}=2$
$x_{1}=\mbox{matrix}(\mbox{rnorm}(n_{1}*p_{1}),n_{1},p_{1})$ | $x_{2}=\mbox{matrix}(\mbox{rnorm}(n_{2}*p_{2},\mbox{mean}=1),n_{2},p_{2})$ | $x_{3}=\mbox{matrix}(\mbox{rnorm}(n_{3}*p_{3}),n_{3},p_{3})$
$\beta_{1}=\mbox{c}(\mbox{rnorm}(s_{1}),\mbox{rep}(0,p_{1}-s_{1}))$ | $\beta_{2}=\mbox{c}(\mbox{rnorm}(s_{2}),\mbox{rep}(0,p_{2}-s_{2}))$ | $\beta_{3}=\mbox{c}(\mbox{rnorm}(s_{3},\mbox{mean}=1),\mbox{rep}(0,p_{3}-s_{3}))$
$y_{1}=x_{1}*\beta_{1}+\mbox{rnorm}(n_{1})$ | $y_{2}=x_{2}*\beta_{2}+\mbox{rnorm}(n_{2})$ | $y_{3}=x_{3}*\beta_{3}+\mbox{rnorm}(n_{3})$
$n_{01}=5$ | $n_{02}=5$ | $n_{03}=5$
$x_{01}=\mbox{matrix}(\mbox{rnorm}(n_{01}*p_{1}),n_{01},p_{1})$ | $x_{02}=\mbox{matrix}(\mbox{rnorm}(n_{02}*p_{2},\mbox{mean}=1),n_{02},p_{2})$ | $x_{03}=\mbox{matrix}(\mbox{rnorm}(n_{03}*p_{3}),n_{03},p_{3})$
$y_{01}=x_{01}*\beta_{1}+\mbox{rnorm}(n_{01})$ | $y_{02}=x_{02}*\beta_{2}+\mbox{rnorm}(n_{02})$ | $y_{03}=x_{03}*\beta_{3}+\mbox{rnorm}(n_{03})$
$D_{1}=\mbox{as.data.frame}(\mbox{cbind}(y_{1},x_{1}))$ | $D_{2}=\mbox{as.data.frame}(\mbox{cbind}(y_{2},x_{2}))$ | $D_{3}=\mbox{as.data.frame}(\mbox{cbind}(y_{3},x_{3}))$
$x_{1}=\mbox{as.data.frame}(x_{1})$ | $x_{2}=\mbox{as.data.frame}(x_{2})$ | $x_{3}=\mbox{as.data.frame}(x_{3})$
Table 2: Simulation datasets with smaller $n$ and large $p$ (using functions from the R language) Cosine | Confomal | Split | Jackknife
---|---|---|---
variable | $1$ | $2$ | $3$ | $1$ | $2$ | $3$ | $1$ | $2$ | $3$
y0 | $0.96$ | $2.05$ | $0.47$ | $0.96$ | $2.05$ | $0.47$ | $0.96$ | $2.05$ | $0.47$
pred | $0.96$ | $2.59$ | $0.95$ | $0.96$ | $2.79$ | $0.9$ | $2.06$ | $0.71$ | $0.98$
predr | $0.88$ | $2.43$ | $0.84$ | $0.79$ | $2.81$ | $0.75$ | $0.67$ | $2.19$ | $0.84$
predrs | $0.88$ | $2.41$ | $0.84$ | $0.91$ | $2.63$ | $0.9$ | $0.69$ | $2.2$ | $0.8$
predl | $0.96$ | $2.59$ | $0.95$ | $0.96$ | $2.59$ | $0.95$ | $2.05$ | $0.71$ | $0.99$
predlr | $0.88$ | $2.43$ | $0.84$ | $0.79$ | $2.81$ | $0.75$ | $0.64$ | $2.08$ | $0.78$
predlrs | $0.9$ | $2.47$ | $0.86$ | $0.91$ | $2.63$ | $0.9$ | $0.65$ | $2.19$ | $0.8$
lo | $-0.36$ | $1.36$ | $-0.36$ | $-0.83$ | $1.01$ | $-0.89$ | $0.22$ | $-1.13$ | $-0.85$
lor | $0.65$ | $0.79$ | $1.08$ | $-1.28$ | $0.36$ | $-0.88$ | $-0.73$ | $0.59$ | $-0.48$
lors | $-0.36$ | $1.08$ | $-0.36$ | $-0.59$ | $0.49$ | $-0.08$ | $-0.63$ | $0.71$ | $-0.42$
lol | $-0.36$ | $1.36$ | $-0.36$ | $-0.81$ | $0.83$ | $-0.81$ | $0.22$ | $-1.12$ | $-0.84$
lolr | $0.5$ | $0.79$ | $1.08$ | $-1.28$ | $0.37$ | $-0.88$ | $-0.73$ | $-0.41$ | $-0.62$
lolrs | $-0.36$ | $1.08$ | $-0.36$ | $-0.59$ | $0.49$ | $-0.08$ | $-0.67$ | $0.71$ | $-0.42$
up | $2.22$ | $3.8$ | $2.22$ | $2.74$ | $4.58$ | $2.69$ | $3.89$ | $2.54$ | $2.82$
upr | $3.08$ | $3.51$ | $3.51$ | $2.96$ | $4.52$ | $2.39$ | $2.08$ | $3.79$ | $2.15$
uprs | $2.08$ | $3.8$ | $2.08$ | $2.39$ | $4.51$ | $1.88$ | $2$ | $3.69$ | $2.03$
upl | $2.22$ | $3.8$ | $2.22$ | $2.72$ | $4.36$ | $2.71$ | $3.88$ | $2.54$ | $2.82$
uplr | $3.08$ | $3.51$ | $3.51$ | $2.96$ | $4.53$ | $2.39$ | $2$ | $4.57$ | $2.19$
uplrs | $2.22$ | $3.8$ | $2.08$ | $2.39$ | $4.51$ | $1.88$ | $1.95$ | $3.57$ | $2.01$
Table 3: Cosine Relevance Note: pred, lo and up denote the prediction, lower limit and upper limit of the interval, respectively. In turn, r correspond to the case when only relevant observations are used while rs appears in cases in which both relevant and simulated observations are included. Finally, l and k refer to the cases in which LASSO regressions and kernel estimations are used, respectively. Percentile | Confomal | Split | Jackknife
---|---|---|---
variable | $1$ | $2$ | $3$ | $1$ | $2$ | $3$ | $1$ | $2$ | $3$
y0 | $0.96$ | $2.05$ | $0.47$ | $0.96$ | $2.05$ | $0.47$ | $0.96$ | $2.05$ | $0.47$
pred | $0.96$ | $2.59$ | $0.95$ | $0.96$ | $2.79$ | $0.9$ | $2.06$ | $0.71$ | $0.98$
predr | $0.84$ | $2.57$ | $0.69$ | $0.43$ | $2.05$ | $0.58$ | $0.67$ | $2.19$ | $0.84$
predrs | $0.86$ | $2.56$ | $0.68$ | $0.7$ | $2.97$ | $0.72$ | $0.69$ | $2.2$ | $0.8$
predl | $0.96$ | $2.59$ | $0.95$ | $0.96$ | $2.59$ | $0.95$ | $2.05$ | $0.71$ | $0.99$
predlr | $0.84$ | $2.57$ | $0.69$ | $0.43$ | $2.02$ | $0.58$ | $0.64$ | $2.08$ | $0.78$
predlrs | $0.87$ | $2.64$ | $0.7$ | $0.7$ | $2.97$ | $0.72$ | $0.65$ | $2.19$ | $0.8$
lo | $-0.36$ | $1.36$ | $-0.36$ | $-0.34$ | $1.5$ | $-0.39$ | $0.22$ | $-1.13$ | $-0.85$
lor | $0.76$ | $0.97$ | $-0.28$ | $-0.21$ | $1.25$ | $-0.59$ | $-0.73$ | $0.59$ | $-0.48$
lors | $-0.25$ | $1.2$ | $-0.41$ | $-0.34$ | $0.88$ | $-0.79$ | $-0.63$ | $0.71$ | $-0.42$
lol | $-0.36$ | $1.36$ | $-0.36$ | $-0.31$ | $1.33$ | $-0.32$ | $0.22$ | $-1.12$ | $-0.84$
lolr | $1.05$ | $1.2$ | $-0.47$ | $-0.26$ | $1.25$ | $-0.6$ | $-0.73$ | $-0.41$ | $-0.62$
lolrs | $-0.18$ | $1.31$ | $-0.41$ | $-0.34$ | $0.88$ | $-0.79$ | $-0.67$ | $0.71$ | $-0.42$
up | $2.22$ | $3.8$ | $2.22$ | $2.25$ | $4.09$ | $2.19$ | $3.89$ | $2.54$ | $2.82$
upr | $2.92$ | $3.49$ | $1.91$ | $2.1$ | $4.09$ | $1.75$ | $2.08$ | $3.79$ | $2.15$
uprs | $1.91$ | $3.94$ | $1.79$ | $2$ | $3.9$ | $2.24$ | $2$ | $3.69$ | $2.03$
upl | $2.22$ | $3.8$ | $2.22$ | $2.22$ | $3.86$ | $2.21$ | $3.88$ | $2.54$ | $2.82$
uplr | $3.29$ | $3.83$ | $1.72$ | $2.17$ | $4.09$ | $1.76$ | $2$ | $4.57$ | $2.19$
uplrs | $1.91$ | $3.94$ | $1.79$ | $2$ | $3.9$ | $2.24$ | $1.95$ | $3.57$ | $2.01$
Table 4: Percentile relevance Note: Here we use the same abbreviations as in
table 3 Short data
---
Percentile | General | Conformal | Split | Jackknife
diffpred | 0.57 | 0.34 | 0.39 | 0.98
diffpredr | 0.25 | 0.28 | 0.21 | 0.26
diffpredrs | 0.33 | 0.27 | 0.48 | 0.25
diffpredl | 0.55 | 0.34 | 0.34 | 0.98
diffpredlr | 0.24 | 0.28 | 0.22 | 0.22
diffpredlrs | 0.35 | 0.31 | 0.48 | 0.26
diffpredk | 0.32 | 0.3 | 0.31 | 0.34
diffpredkr | 0.24 | 0.2 | 0.24 | 0.28
diffpredkrs | 0.25 | 0.19 | 0.27 | 0.28
%pred | 60.25 | 42.49 | 42.38 | 95.87
%predr | 30.52 | 27.57 | 26.09 | 37.92
%predrs | 34.42 | 26.03 | 41.81 | 35.41
%predl | 60.31 | 42.49 | 42.56 | 95.89
%predlr | 29.28 | 27.57 | 26.65 | 33.62
%predlrs | 35.9 | 29.26 | 41.91 | 36.54
%predk | 38.37 | 36.05 | 36.49 | 42.58
%predkr | 26.22 | 23.44 | 27.83 | 27.39
%predkrs | 27.79 | 23.29 | 32.62 | 27.45
int | 2.93 | 2.53 | 2.59 | 3.67
intr | 2.56 | 2.29 | 2.49 | 2.89
intrs | 2.62 | 2.37 | 2.8 | 2.69
intl | 2.91 | 2.53 | 2.53 | 3.66
intlr | 2.8 | 2.35 | 2.54 | 3.5
intlrs | 2.58 | 2.31 | 2.8 | 2.64
intk | 2.56 | 2.53 | 2.49 | 2.65
intkr | 2.95 | 2.56 | 2.61 | 3.66
intkrs | 2.45 | 2.29 | 2.65 | 2.41
ab | 0.6 | 0.63 | 0.65 | 0.52
abr | 0.61 | 0.71 | 0.59 | 0.53
abrs | 0.55 | 0.58 | 0.55 | 0.53
abl | 0.6 | 0.63 | 0.63 | 0.52
ablr | 0.62 | 0.76 | 0.59 | 0.5
ablrs | 0.55 | 0.59 | 0.55 | 0.51
abk | 0.62 | 0.61 | 0.62 | 0.63
abkr | 0.57 | 0.63 | 0.54 | 0.55
abkrs | 0.54 | 0.54 | 0.51 | 0.58
Table 5: Short data percentile Note: diffpred denotes the differences between
the predicted and the true values. %pred indicates the same difference but
relative to the true value. int is the length of the interval (the upper minus
the lower limits). Finally, ab represents the size of the error relative to
the length of the interval. Short data
---
Cosine | General | Conformal | Split | Jackknife
diffpred | 0.57 | 0.34 | 0.39 | 0.98
diffpredr | 0.31 | 0.28 | 0.4 | 0.26
diffpredrs | 0.29 | 0.27 | 0.35 | 0.25
diffpredl | 0.55 | 0.34 | 0.34 | 0.98
diffpredlr | 0.3 | 0.28 | 0.4 | 0.22
diffpredlrs | 0.3 | 0.29 | 0.35 | 0.26
diffpredk | 0.32 | 0.3 | 0.34 | 0.31
diffpredkr | 0.6 | 0.8 | 0.57 | 0.44
diffpredkrs | 0.93 | 1.1 | 0.58 | 1.11
%pred | 60.25 | 42.49 | 42.38 | 95.87
%predr | 37.08 | 35.18 | 38.15 | 37.92
%predrs | 37.27 | 34.79 | 41.6 | 35.41
%predl | 60.31 | 42.49 | 42.56 | 95.89
%predlr | 35.65 | 35.18 | 38.15 | 33.62
%predlrs | 38.18 | 36.4 | 41.6 | 36.54
%predk | 38.41 | 36.05 | 42.7 | 36.49
%predkr | 89.13 | 102.14 | 99.06 | 66.2
%predkrs | 144.67 | 153.04 | 99.05 | 181.92
int | 3.26 | 2.53 | 3.58 | 3.67
intr | 3.1 | 2.53 | 3.89 | 2.89
intrs | 2.74 | 2.53 | 2.99 | 2.69
intl | 3.24 | 2.53 | 3.52 | 3.66
intlr | 3.33 | 2.58 | 3.89 | 3.5
intlrs | 2.74 | 2.58 | 2.99 | 2.64
intk | 2.56 | 2.53 | 2.64 | 2.49
intkr | 2.93 | 2.86 | 3.41 | 2.52
intkrs | 2.91 | 3.37 | 2.68 | 2.67
ab | 0.59 | 0.63 | 0.61 | 0.52
abr | 0.65 | 0.88 | 0.55 | 0.53
abrs | 0.57 | 0.59 | 0.6 | 0.53
abl | 0.58 | 0.63 | 0.6 | 0.52
ablr | 0.64 | 0.87 | 0.55 | 0.5
ablrs | 0.57 | 0.6 | 0.6 | 0.51
abk | 0.62 | 0.61 | 0.63 | 0.62
abkr | 0.62 | 0.37 | 0.54 | 0.94
abkrs | 0.73 | 0.62 | 0.58 | 0.99
Table 6: Short data cosine Note: Here we use the same abbreviations as in
table 5 Long data
---
Cosine | General | Conformal | Split | Jackknife
diffpred | 2.8 | 2.47 | 2.58 | 3.34
diffpredr | 2.75 | 2.27 | 2.64 | 3.34
diffpredrs | 2.8 | 2.29 | 2.76 | 3.34
diffpredl | 2.81 | 2.48 | 2.6 | 3.35
diffpredlr | 2.73 | 2.29 | 2.61 | 3.28
diffpredlrs | 2.78 | 2.34 | 2.69 | 3.33
diffpredk | 2.05 | 1.5 | 2.32 | 2.33
diffpredkr | 1.91 | 1.7 | 1.87 | 2.17
diffpredkrs | 2.1 | 1.77 | 2.39 | 2.14
%pred | 179.03 | 165.13 | 105.96 | 266.02
%predr | 168.86 | 157.93 | 115.92 | 232.74
%predrs | 183.81 | 157.16 | 165.13 | 229.12
%predl | 176.34 | 165.17 | 98.11 | 265.75
%predlr | 173.31 | 135.75 | 109.4 | 274.78
%predlrs | 174.9 | 160.28 | 136.48 | 227.92
%predk | 161.98 | 132.38 | 183.64 | 169.93
%predkr | 131.49 | 166.7 | 109.41 | 118.35
%predkrs | 154.95 | 126.67 | 215.02 | 123.15
int | 5.13 | 5.86 | 5.96 | 3.57
intr | 5.03 | 5.51 | 5.95 | 3.64
intrs | 4.31 | 5.87 | 3.52 | 3.54
intl | 5.14 | 5.69 | 6.38 | 3.36
intlr | 6.78 | 5.76 | 6.24 | 8.33
intlrs | 4.31 | 5.94 | 3.52 | 3.48
intk | 5.43 | 4.12 | 6.71 | 5.45
intkr | 6.05 | 5.27 | 5.58 | 7.31
intkrs | 5.97 | 1.77 | 9.41 | 6.73
ab | 0.6 | 0.42 | 0.43 | 0.94
abr | 0.59 | 0.41 | 0.45 | 0.92
abrs | 0.71 | 0.39 | 0.8 | 0.94
abl | 0.61 | 0.44 | 0.41 | 1
ablr | 0.41 | 0.4 | 0.42 | 0.4
ablrs | 0.75 | 0.39 | 0.78 | 1.09
abk | 0.38 | 0.36 | 0.35 | 0.43
abkr | 0.33 | 0.33 | 0.35 | 0.3
abkrs | 0.63 | 1.21 | 0.36 | 0.32
Table 7: Long data cosine Note: Here we use the same abbreviations as in table
5