# Probing out-of-distribution generalization in machine learning for materials

Kangming Li  0000-0003-4471-8527<EMAIL_ADDRESS>Department of
Materials Science and Engineering, University of Toronto, 27 King’s College
Cir, Toronto, ON, Canada.  Andre Niyongabo Rubungo  0000-0003-3608-2039
Vertaix, Department of Computer Science, Princeton University, Princeton, NJ,
08544, USA.  Xiangyun Lei  Toyota Research Institute, 4440 El Camino Real, Los
Altos, California 94022, USA.  Daniel Persaud  Department of Materials Science
and Engineering, University of Toronto, 27 King’s College Cir, Toronto, ON,
Canada.  Kamal Choudhary  0000-0001-9737-8074 Material Measurement Laboratory,
National Institute of Standards and Technology, 100 Bureau Dr, Gaithersburg,
MD, USA.  Brian DeCost  0000-0002-3459-5888 Material Measurement Laboratory,
National Institute of Standards and Technology, 100 Bureau Dr, Gaithersburg,
MD, USA.  Adji Bousso Dieng  0000-0001-5687-3554 Vertaix, Department of
Computer Science, Princeton University, Princeton, NJ, 08544, USA.  Jason
Hattrick-Simpers  0000-0003-2937-3188<EMAIL_ADDRESS>Department of Materials Science and Engineering, University of Toronto, 27
King’s College Cir, Toronto, ON, Canada.  Acceleration Consortium, University
of Toronto, 27 King’s College Cir, Toronto, ON, Canada.  Vector Institute for
Artificial Intelligence, 661 University Ave, Toronto, ON, Canada.  Schwartz
Reisman Institute for Technology and Society, 101 College St, Toronto, ON,
Canada.

###### Abstract

Scientific machine learning (ML) endeavors to develop generalizable models
with broad applicability. However, the assessment of generalizability is often
based on heuristics. Here, we demonstrate in the materials science setting
that heuristics based evaluations lead to substantially biased conclusions of
ML generalizability and benefits of neural scaling. We evaluate generalization
performance in over 700 out-of-distribution tasks that features new chemistry
or structural symmetry not present in the training data. Surprisingly, good
performance is found in most tasks and across various ML models including
simple boosted trees. Analysis of the materials representation space reveals
that most tasks contain test data that lie in regions well covered by training
data, while poorly-performing tasks contain mainly test data outside the
training domain. For the latter case, increasing training set size or training
time has marginal or even adverse effects on the generalization performance,
contrary to what the neural scaling paradigm assumes. Our findings show that
most heuristically-defined out-of-distribution tests are not genuinely
difficult and evaluate only the ability to interpolate. Evaluating on such
tasks rather than the truly challenging ones can lead to an overestimation of
generalizability and benefits of scaling.

##  I Introduction

Machine learning (ML) has emerged as an important tool in accelerating
scientific discovery [1, 2, 3, 4]. This transition to data-driven science is
epitomized by the development of scientific ML that seeks to build
generalizable models capable of broad applicability [5]. In chemical and
materials sciences, recent studies have been focused on developing universal
or foundational deep learning models, which are suggested to achieve
unprecedented levels of out-of-distribution (OOD) generalizations towards
unseen materials that are dissimilar to the training data [6, 7, 8, 9, 10,
11].

However, a critical issue that has been overlooked is the potential biases in
selecting OOD tasks to demonstrate generalizability. OOD tasks are often
defined based on simple heuristics, which, due to their subjective nature,
vary between studies and even lead to contradicting interpretation of
generalizability. For instance, the generalization to structures with 5+
elements, despite their omission from training, was used to showcase the
emergent capability of deep learning models [6]. However, this perspective has
been contested with the argument that such generalizations are anticipated
from the physical heuristic that interactions in higher-order systems can be
inferred from lower-order ones [12, 13, 14, 15]. This discrepancy underscores
a broader lack of agreement and discussion on what constitutes a genuinely
challenging OOD task. Indeed, if the OOD test set falls within the training
domain, conclusions in the superiority of a state-of-the-art ML architecture
and the real benefits of model scaling may pertain only to interpolation
capabilities rather than true extrapolation.

In this work, we conduct a systematic examination of the performance of
various ML methods across over 700 OOD tasks within large materials datasets.
These tasks are specifically designed to challenge common heuristics based on
chemistry or structural symmetry. Our findings indicate that existing models,
including those as simple as tree ensembles, exhibit robust generalization
across most tasks that feature new chemical or structural groups absent from
the training data. Our results highlight that most of these heuristics-based
criteria do not constitute truly challenging tasks for ML models. By analyzing
the representation spaces of materials, we reveal that test data from well
performing tasks largely reside within the training domain, whereas those from
poorly performed tasks do not. Furthermore, we find that OOD generalization
performance does not necessarily follow traditional scaling laws [16, 17].
Notably, scaling up training set size or training time leads to marginal
improvement or even degradation in the generalization performance for those
challenging OOD tasks. These findings indicate that the purported benefits of
scaling and emergent generalizability could be considerably overstated, due to
domain misidentification driven by human bias that confounds regimes of
interpolation and extrapolation.

##  II Results

###  II.1 Evaluation setup

ML models are usually evaluated for their in-distribution (ID) performance by
using random train-test split of the whole dataset. In this work, we assess
model performance on OOD tasks, where a specific group of materials serves as
the test set while the remainder of the dataset forms the training data. We
primarily consider six criteria for defining OOD test data: (1) materials
containing element $X$, (2) materials containing any element in the period
$X$, (3) materials containing any element in the group $X$, (4) materials of
space group $X$, (5) materials of point group $X$, (6) materials of crystal
system $X$. We consider these leave-one-$X$-out tasks for all possible values
of $X$ but exclude those with fewer than 200 test samples.

Three ab initio-derived materials databases have been selected for this
evaluation: Joint Automated Repository for Various Integrated Simulations
(JARVIS) [18, 19], Materials Project (MP) [20], and the Open Quantum Materials
Database (OQMD) [21]. These databases have different data distributions, hence
ensuring a robust and generalized conclusion of OOD performance. The
combination of datasets and grouping criteria leads to over 700 OOD tasks. OOD
performance is evaluated by training a representative set of ML models,
including (1) random forest (RF) [22] and XGBoost (XGB) [23] models with
Matminer descriptors [24]. (2) single neural network with Gaussian multipole
(GMP) expansion on electron density [25] (3) atomistic line graph neural
network (ALIGNN) [26], and (4) the large language model (LLM) based LLM-Prop
with crystal text descriptions [27]. It therefore covers not only common ML
architectures from lightweight tree ensembles and neural networks to graph
neural networks and transformer-based large language models, but also distinct
input representations from human-devised descriptors and force-field-like
features to crystal graphs and text.

###  II.2 Chemistry-based OOD generalization

In the following, our analysis primarily focuses on formation energy data,
given its fundamental importance in materials science. Table 1 shows the in
distribution (ID) performance obtained from a random 8:2 train-test split of
the whole dataset, providing a baseline for comparison with OOD performance.
We utilize two complementary performance metrics: mean absolute error (MAE)
and coefficient of determination ($R^{2}$). While MAE measures the expected
error on the original physical scale, its interpretation can be scale-
dependent and less intuitive for assessing the goodness of predictions. In
contrast, $R^{2}$ is a dimensionless accuracy measure with a value range
between 1 (for a perfect model) and infinitely negative (for arbitrarily bad
models) and can be conveniently compared across different OOD test sets.

Table 1: In-distribution performance for the formation energy prediction. Models are arranged in the ascending order of MAEs of the MP dataset from left to right. Best performance is highlighted in bold. The MP, JARVIS, and OQMD datasets contain 146k, 76k, and 1M entries, respectively.  Metric | Dataset | ALIGNN | GMP | LLM-Prop | XGB | RF
---|---|---|---|---|---|---
MAE (eV/at) | MP | 0.033 | 0.052 | 0.063 | 0.078 | 0.090
JARVIS | 0.036 | 0.081 | 0.068 | 0.074 | 0.099
OQMD | 0.020 | 0.038 | 0.045 | 0.070 | 0.065
$R^{2}$ | MP | 0.996 | 0.992 | 0.981 | 0.979 | 0.970
JARVIS | 0.995 | 0.985 | 0.982 | 0.981 | 0.968
OQMD | 0.998 | 0.995 | 0.995 | 0.987 | 0.985
Figure 1: Leave-one-element-out performance on Materials Project. (a) Leave-
one-element-out performance across the periodic table for the XGB and ALIGNN
models. Despite the presence of negative $R^{2}$ for H and F, the colorbar
range is bounded between 0 and 1 for better visibility. The OOD tests are
performed only for elements with more than 200 data. Results for other models
and datasets are provided in Supplemental Material. The periodic tables are
plotted using pymatviz [28]. (b) Violin plots of SHAP contributions from
compositional and structural features for selected left-out elements. (c)
Number of leave-one-element-out tasks with $R^{2}$ higher than a given
threshold as a function of threshold. LLM-Prop is not included due to the
large number of tasks and high training cost. (d) Comparison of performance
between ML models for selected leave-one-element-out tasks.

Figure 1(a) illustrates the leave-one-element-out generalization performance
of the ALIGNN and XGB models on the MP dataset, demonstrating robust OOD
performance across much of the periodic table. Notably, 85 % of the tasks for
the ALIGNN model and 68 % for the simpler XGB model achieved $R^{2}$ scores
above 0.95. This broad generalization capability is surprising, given that the
training set does not contain any information regarding the bonding between
the left-out element and other atoms. Our results suggest that effective OOD
generalization across chemistry may be more achievable than previously
assumed, for both low- and high-capacity models.

Tasks with low $R^{2}$ scores are mainly associated with nonmetals such as H,
F, O. Nonetheless, ML models may still prove useful in ranking materials even
in the worst performing tasks. As shown in the parity plots, the ALIGNN model
systematically overestimates the formation energies of H compounds (or O
compounds), with a Pearson correlation coefficient of 0.7 (or 0.9). The poor
$R^{2}$ is therefore a result of systematic biases in the OOD predictions.
Such bias is also shared among other ML models: all ML models considered here
overestimate formation energies of H and F compounds; the three deep learning
models tend to overestimate formation energies of O compounds, while tree
models tend to overestimate low-energy and underestimate high-energy O
compounds. Similar systematic biases, which can be addressed by simple linear
corrections, have also been reported elsewhere for OOD property predictions
[29, 30]. Future investigation to understand and mitigate these biases can be
useful for developing chemically transferable models.

The causes behind the poor OOD performance for nonmetals, particularly H, F,
and O, have been examined. Initially, we considered training set size as a
potential factor, but this was ruled out based on the weak correlation between
training set size and OOD performance, evidenced by a Spearman’s rank
correlation coefficient of -0.2. Despite investigating various element
attributes, none could fully account for the $R^{2}$ trends. For example,
while elements with low $R^{2}$ scores tend to be more electronegative or
positioned at the boundary of the periodic table, exceptions such as S, Cl, Br
(which are electronegative) and Cs, Bi (positioned near the corner of the
periodic table) exhibit high ALIGNN $R^{2}$ scores above 0.96, suggesting
other influencing factors.

The biases leading to this poor OOD performance could be of either
compositional or structural origin. The former is related to chemical
dissimilarity between elements, whereas the latter refers to the failure mode
where the test set associated with a given element contains materials that are
structurally distinct from the training set. We propose a SHAP-based [31]
(SHapley Additive exPlanations) method to identify the sources of biases. This
involves training a second model to correct the model for each leave-one-
element-out task and evaluating the contributions from compositional and
structural features to the corrections (see Method). Fig. 1(b) shows the
violin plots of the feature contributions from the XGB model for the test data
in selected leave-one-element-out tasks. Compared to the well-performing case
of Cl, the compositional contributions in the cases of H, F, and O are much
stronger in magnitude than structural contributions. The predominately
negative compositional contributions are consistent with the corrections
needed for the overestimation shown in the parity plots. This analysis
therefore reveals that the poor OOD performance is mainly of chemical origin,
namely the chemistry-property relationship is drastically different and cannot
be readily transferred from other chemical systems.

Figure 2: Leave-one-period-out and leave-one-group-out performance on
Materials Project. The inset in (a) shows the parity plot for the ALIGNN
predictions for the period 5 (axis range: -5 eV/atom to 5 eV/atom).

To compare the overall performance between different models, we count for each
model the number of leave-one-element-out tasks with an accuracy better than a
given threshold. As shown in Fig. 1(c), models with better ID performance have
a better overall performance in the leave-one-element-out tasks. However, this
correlation does not imply uniform superiority of any model across all tasks.
Figure 1(d) details performance comparisons in 10 representative tasks,
highlighting nuanced outcomes. In the most challenging tasks involving H, F,
and O, the RF model, despite its lower ID performance, ranks as the best or
second best. Conversely, for moderately challenging tasks (P, I, and Li), the
GMP model outperforms others. Meanwhile, in less demanding tasks (Bi, Fe, and
Rb) where MAEs are below 0.1 eV/atom, the ALIGNN model achieves a much lower
MAE than the other models.

Interestingly, the LLM-Prop model, despite its good ID performance, is the
worst performing model in most of the leave-one-element-out tasks. A possible
cause might be the lack of inductive bias related to the element similarity,
which is encoded in the representation schemes of other models. For instance,
ALIGNN and tree ensembles represent atoms by basic attributes such as
electronegativity, period, and group. In this way, the models can properly
recognize the correlation between the left-out element and the other elements
in the training set. By contrast, language models treat each element as a text
token in a manner equivalent to one-hot encoding. Element similarity may not
be encoded in the initial text embedding and is not guaranteed to be learned
for the left-out element during the training. A potential direction to improve
chemical transferability of LLMs is therefore to include additional chemical
information in text descriptions during pretraining or fine-tuning.
Alternatively, language models specialized for materials and chemistry
applications might explore custom tokenization schemes for elemental symbols
and chemical formulae that directly incorporate aspects of chemical
similarity.

It is worth noting that many representation schemes that use one-hot encoding
of elements or element-specific functions [32] are in principle incapable of
any generalization to new chemistry. By contrast, the good generalization
across chemistry shown in Fig. 1 highlights the importance of describing
elements with a transferable representation. In this way, an element and its
compounds can be seen as lying within the interpolation region in the chemical
space of the training data in the leave-one-element-out tasks.

Appropriate representations can achieve chemical transferability even for
small training sets with limited chemical diversity. To prove this point, we
create a dataset from OQMD by restricting the chemistry to six elements from
Cr to Cu. The leave-one-element-out $R^{2}$ scores of the XGB model within
this dataset (about 1200 structures) range between 0.87 (for Cu) to 0.93 (for
Fe), comparable to the in-distribution $R^{2}$ score of 0.93, demonstrating
that generalization across chemistry does not necessarily require a big model
or large dataset.

Fig. 2(a) shows the more challenging OOD tasks where an entire period of
elements is excluded from the training set. The worst OOD performance is found
for the first period, which contains H, and the second period, which contains
many nonmetals with poor leave-one-element-out performance. Overall, the OOD
performance tends to be better for the periods in the middle of the periodic
table, where ML models can learn from the adjacent periods. For the OOD tasks
of the fourth to sixth periods, the $R^{2}$ of the ALIGNN model can achieve a
$R^{2}$ score of 0.95 to 0.97, which is surprisingly good given the number of
elements excluded from the training set. Interestingly, the XGB and RF models
perform similarly or even better than deep learning models in most of the
tasks, highlighting the robustness of tree ensembles in challenging
generalization tasks.

Similar conclusions are also found in the leave-one-group-out performance
shown in Fig. 2(b). First, the OOD performance also tends to be better for the
groups in the middle of the periodic table. Despite its superior leave-one-
element-out performance for elements in groups 3 to 5, the ALIGNN model is
outperformed by the tree ensembles in these leave-one-group-out tasks. By
contrast, the GMP model demonstrates the best or close-to-best performance in
these tasks. The difference in the model performance ranking in the leave-one-
element/period/group-out tasks indicates that these ML models have distinct
advantages in different types of OOD tasks. Future investigation in
understanding their respective advantages can be beneficial for designing
models with better generalizability.

In the JARVIS and OQMD datasets, we also find good generalization performance
across the majority of the periodic table and poor performance for a few cases
of nonmetal elements. Therefore, the trend shown in Fig. 1 and 2 is not a
unique artifact of the specific biases in data distribution or computational
settings of the MP dataset, nor of the training set size which will be
discussed later. Instead, it reveals that the existing featurization schemes
can enable transferability between most of the chemical systems, though they
may still fail to capture the essential characteristics for certain nonmetals,
which leaves rooms for future improvement in materials featurization.

###  II.3 Structure-based OOD generalization

Figure 3:  Structure-based OOD performance on Materials Project. (a) Number of
leave one space group out tasks with $R^{2}$ higher than a given threshold as
a function of threshold. The horizontal dashed line indicates the total number
of tasks. LLM-Prop is not examined due to the large number of tasks and high
training cost. (b) Leave one crystal system out performance in the ascending
order of symmetry from left to right. (c) Leave one prototype-like out
performance. Structures with fingerprints similar to a prototype are grouped
together. (d) Leave one point group out performance in the ascending order of
symmetry from left to right.

Structures are usually grouped by their crystallographic symmetries. Here we
consider three basic symmetry attributes for structural grouping: space
groups, point groups, and crystal systems, with each being a subcategory of
the next. The structure-based OOD performance on the MP formation energy
dataset are shown in Fig. 3.

The OOD performance of the structure-based tasks is much better than that of
the chemistry-based OOD tasks. Fig. 3(a) shows the overall OOD performance on
all of the leave one space group out tasks. As the best performing model,
ALIGNN achieves an $R^{2}$ of above 0.95 in all of the tasks with 88 % of them
having an $R^{2}$ score above 0.98. The GMP model is the second best
performing one followed by the XGB model. The LLM-Prop model is not examined
here due to the high training cost required for a large number of tasks, but
it is included for the leave one point group out and leave one crystal system
out tasks in Fig. 3(b) and (c). All models demonstrate a good prediction
accuracy with an $R^{2}$ of above 0.95 (or a MAE of below 0.15 eV/atom) in
most of these tasks. The only cases with relatively poor performance are
associated with the triclinic crystal system (and its subcategory point groups
1 and $\overline{1}$), which has minimal symmetry constraint and is
structurally more diverse, posing greater challenge when generalized from
structures with higher symmetry.

Our results suggest that symmetry attributes, which are often used to explain
generalization performance or devise active learning algorithms [33, 34, 35,
9, 36], may not be effective in finding data truly dissimilar from the
training set. Indeed, two symmetrically different crystals can be structurally
very similar. For example, a slight distortion of a cubic lattice can
transform the structure into another lattice, as in the famous cubic-to-
tetragonal Martensitic transformation in steels. Therefore, the training set
still contains similar structures even if all cubic structures are removed. To
mitigate this information leakage and make the generalization task more
challenging, we create new test sets by using the fuzzy prototype-matching
algorithm implemented in the robocrystallographer package [37] to group
prototype-like structures. Fig. 3(d) shows the OOD performance on a few leave
one prototype-like out tasks. Interestingly, the OOD performance in these
tasks is still good, especially for deep learning models. The ALIGNN and GMP
models achieve $R^{2}$ scores above 0.95 in all tasks except for Heusler-like
structures, for which the ALIGNN model can still achieve an $R^{2}$ of 0.86
and an MAE of 0.08 eV/atom. The XGB model still perform reasonably well in
these tasks, though its MAEs are 100 % to 200 % higher than those of the
ALIGNN model.

Overall, the ranking of model performance in structure-based OOD tasks closely
follows the ID performance ranking, in contrast to the situation in chemistry-
based OOD tasks. In particular, the ALIGNN model systematically outperforms
other models in all the structure-based tasks, highlighting the effectiveness
of graph neural networks in capturing complex structural patterns. In
addition, the LLM-Prop model also exhibits much better generalization across
structural groups than across chemical groups.

###  II.4 Inspecting materials representations

Figure 4: UMAP two-dimensional projection of materials representations for
domain identification. (a-b) UMAP plots of ALIGNN embeddings learned from the
leave-Mg-out and leave-O-out tasks in JARVIS. (c) UMAP plot of ALIGNN
embeddings learned from the leave-H-out task in OQMD. (d) UMAP plot of ALIGNN
embeddings learned by leaving out structures with 5 or more elements in MP.
(e) UMAP plot of the XGB descriptors for the leave-period-5-out task in MP.
(f) Absolute errors (left Y axis) of test data as functions of kernel-density
estimates of training data for the UMAP plot of (e); the solid line denotes
the MAEs (right Y axis) for different density intervals. In all UMAP plots,
the training data, in-domain test data, and out-of-domain test data are marked
in grey, blue, and red, respectively; clusters of out-of-domain test data are
circled out in (a-c); the $R^{2}$ scores are indicated for the in-domain, out-
of-domain, and all-domain test data.

The systematic investigation of various OOD tasks reveals that in most of the
cases, generalization beyond the distribution of the training data is not
challenging. It therefore seems to suggest that ML models are good at
extrapolating across chemical and structural groups, in contrast to the recent
findings that ML models generalize poorly even between different versions of
the same databases [33, 38]. To reconcile this apparent contradiction, we
propose to inspect materials not from the perspective of a single human-
defined attribute (e.g., element identity), but from the high-dimensional
representation space of materials. Furthermore, we propose to define and
distinguish two types of OOD. A test set is considered statistically OOD if
the statistical distribution of one or more attributes differs from the
training set, whereas a test set is considered representationally OOD if the
test data are outside the region of the training set when viewed from the
high-dimensional representation space. With this definition, all the OOD tasks
discussed previously are statistically OOD (which will be simply referred to
as OOD) but not necessarily representationally OOD.

Uniform Manifold Approximation and Projection (UMAP) [39] is used to project
the high-dimensional materials representations to a two-dimensional plane. For
deep learning models, materials representations (embeddings) are automatically
learned during the training; embeddings of new materials can then be derived
by using the trained models. For each OOD task, we apply a Gaussian kernel to
the training data and evaluate the kernel density estimate of the training
data for every test data point in the UMAP embedding space. Test data points
with high/low density estimates are considered as representationally ID/OOD
data.

Figure 5: Effects of training time and training set size on performance. (a-c)
Training loss, ID test loss, and OOD test loss as functions of training
round/epoch on JARVIS. Models and left-out elements are indicated in the top
of the plots. (d-e) ALIGNN performance on Materials Project and OQMD,
respectively. Solid lines with markers denotes the MAEs as functions of
training set size for 3 representative OOD tasks: leave F out, leave H out,
and leave out structures with $\geq$ 5 elements. Dash line marks the ID MAE as
the baseline for comparison. (f) Distribution of ALIGNN prediction errors on
H-containing structures in the leave-H-out task in OQMD. Results for the
models trained on 1 %, 33 %, and all of the structures without H are shown.
Similar analysis for other tasks and ML models is provided in Supplemental
Material.

Fig. 4(a-b) show the UMAP plots of ALIGNN embeddings learned from the leave-
Mg-out and leave-O-out tasks in the JARVIS dataset. In the former case, 76 %
of Mg-containing structures are well covered by the training data (i.e., with
an embedding space probability density above the threshold), whereas 24 % of
Mg-containing structures are representationally OOD with a much lower $R^{2}$
score. In the latter task, 87 % of O-containing structures are
representationally OOD with an $R^{2}$ close to 0, whereas there are still
some in-domain O-containing structures with an $R^{2}$ above 0.9. Fig. 4(c)
shows another example from the leave-H-out task in the OQMD dataset, where
kernel-density estimates provide again good resolution in separating
representationally ID and OOD test data: 13 % of H-containing structures are
identified as within the training domain with an $R^{2}$ score of 0.8. The
identification of poorly predicted test data in well-performing OOD tasks and
well-predicted test data in poorly-performing OOD tasks demonstrates the
effectiveness of the proposed method for classifying representational domains
of OOD samples.

Representational domain identification is important for interpreting
generalization capabilities. Fig. 4(d) shows the case where the test set
contains structures with 5 or more elements. It has been used as an example to
demonstrate the emergent ability of deep learning models [6]. This perspective
was contested as the same capability was also found in simple boosted trees
[12]. In fact, our UMAP analysis shows that 92 % the test data are identified
as lying within the representational domain of the training data, therefore
this test only reflects the interpolation capacity of ML models. This example
highlights not only the risk of misinterpreting generalization capabilities,
but also the need for designing more challenging tests that can truly
demonstrate state-of-the-art.

Domain identification can be achieved not only with the learned embeddings but
also with the descriptors used in tree ensembles. An example is shown in Fig.
4(e), where the test data are structures that contain any element in the fifth
period. Based on kernel-density estimates, 93 % of the test data are
identified as representationally ID data, thereby explaining the good
generalization performance for this seemingly difficult task.

It is worth noting that density estimates do not perfectly correlate with
prediction errors, as shown in Fig. 4(f). However, from a statistical
perspective there is a higher proportion of test data with high prediction
errors in the low-density region than in the high-density region. Indeed, when
calculating the MAE for test data in different density intervals, we find that
the MAE tends to decrease with increasing density estimate (dotted line). The
density estimate can still enable us to statistically distinguish between the
well-predicted in-domain data and poorly-predicted out-of-domain data from an
out-of-distribution test set. Future work will be needed to improve domain
identification and uncertainty measures with better correlation with OOD
prediction errors.

###  II.5 Learning curves for out-of-distribution generalization

Domain identification in the materials representation space is important in
correctly interpreting the benefits of the neural scaling strategy for
materials discovery applications. This strategy is based on the so-called
neural scaling laws, where test errors can be consistently reduced by
increasing model size, dataset size, and the amount of compute used for
training deep learning models. Recent work on deep learning models for
materials follows this rationale to train foundational or universal neural
networks [6, 10, 40, 8, 9]. However, neural scaling laws have been empirically
tested only for the ID generalization [16, 41]. Here we perform an extensive
examination on various OOD tasks and, intriguingly, find distinct scaling
effects among the OOD tasks. Specifically, neural scaling laws are found to be
valid only for representationally ID generalization, while the scaling effects
are marginally beneficial or even adverse for representationally OOD
generalization.

Fig. 5(a-c) shows the learning curves of MAE versus training time in the
leave-one-element-out tasks for F, O, and H. The training and the ID set
contain 80 % and 20 % of the structures without the left-out element, whereas
the OOD test set contains the structures with the left-out element. For the
leave-F-out task, the training and ID test MAEs decrease with increasing
training time of the XGB model. However, the OOD MAE decreases only during the
initial 20 training rounds, then continuously increases with more training.

A similar pattern is observed in the leave-O-out and leave-H-out tasks. In
fact, for representationally OOD tasks, the test MAE tend to increase or
remain constant with increasing training time beyond around 20 epochs.
However, for OOD tasks that are representationally ID, the OOD MAE has a
tendency similar to the ID MAE, namely the error decreases with increasing
training time. As discussed in Fig. 4, a key distinction between the easy and
difficult OOD tasks is whether the test data lies outside the training domain
in the representation space. Thus, our findings indicate that the response to
scaling up training efforts may vary significantly based on whether the test
data falls inside or outside the training domain.

Effects of training set size are also investigated. Take the leave-F-out task
for example, models are trained with different sizes of the training sets
randomly sampled from 80 % of the structures without F, and evaluate
performance on the held-out ID test set (20 % of the structures without F) and
the OOD test set (structures with F).

Fig. 5(d) shows the ALIGNN learning curves for the ID and OOD performance for
3 representative tasks in Materials Project. As the ID learning curves for
different OOD tasks are found to be almost the same, we show only one learning
curve for the ID performance (dashed line) for simplicity. The ID performance
consistently improves with increasing training set size, which is however not
always the case for OOD performance. For the leave-F-out task, increasing the
training set size from $10^{2}$ to $10^{5}$ only further degrades the already
bad OOD performance. For the leave-H-out task, we find a “V” shaped learning
curve for the performance on H-containing structures, where the MAE first
decreases and then increases with increasing training set size. This violation
of neural scaling laws is also found in other representationally OOD tasks
(Supplemental Material). By contrast, we find a consistent decrease in MAE
with increasing training set size, for the case where the test data are
structures with 5 or more elements. This accordance with neural scaling laws
is also found in other easy OOD test sets that contain mainly
representationally ID data (Supplemental Material).

To avoid establishing our conclusion based on a single dataset, we also
perform the same experiments on the OQMD dataset, which is about an order of
magnitude larger than the Materials Project dataset. As shown in Fig. 5(e), we
find the same trend where increasing training set size only leads to marginal
improvement or even degradation for test data that are outside the training
domain. It should be noted that the marginal improvement is not an artifact of
the limited model capacity, considering the continuous decrease in the ID test
MAE with increasing training set size. In addition, similar learning curves
are also found for other ML models and the JARVIS dataset (Supplemental
Material).

It is particularly striking to observe the 5-fold increase in the OOD MAE of
the leave-H-out task, when the training set size is increased from $10^{4}$ to
$10^{6}$. Fig. 5(f) shows the distributions of the corresponding prediction
errors on H-containing structures. When trained on only 1 % of structures
without H, the error distribution is centered around zero with a relatively
narrow spread. With increasing training set size, however, the errors shift
towards positive values and exhibit a larger spread. This systematic
overestimation is a reflection of biases in models that are enhanced by
training on more data from a given domain. The leave-H-out task can therefore
serve as an example of losing generalizability due to an overfitting to the
training domain, for which future work is needed to identify its root cause
and mitigation solutions.

Our results challenge previous interpretations of generalizability. Figure
5(d-e) include the learning curves for the scenario where the test set
comprises structures with five or more elements, while the training set
includes those with fewer than five. Contrary to prior claims suggesting such
generalization is an emergent capability enabled by neural scaling [6], our
analysis indicates that this may merely reflect in-domain generalization.
According to the literature, a sudden improvement in performance above a
certain threshold is required to classify an ability as emergent [42];
however, this criterion is not met in the learning curve. Additionally, the
OOD MAE being lower than the ID test MAE further supports our findings from
Figure 4(e) that the test data lie well within the representational training
domain. Thus, the observed effective generalization in these cases does not
demonstrate the models’ true generalizability but rather the interpolation
within the representational training domain, where errors decrease as the test
data become increasingly well-covered by the training data.

##  III Discussion

Our extensive OOD examinations show that existing models exhibit remarkable
generalization capabilities to new chemical or structural groups outside the
training set. We explain these observations by inspecting materials
representation, revealing that many OOD test sets are actually within the
representational domain of training data, whereas only a handful of
challenging OOD tasks contain significant amounts of test data outside the
representational training domain. Intriguingly, the performance on these
representationally OOD tasks does not adhere to conventional neural scaling
laws, suggesting that the benefits of scaling for OOD generalization may be
overstated or misinterpreted.

Our findings underscore the predictive power of ML but also reveal that many
OOD tasks commonly considered difficult based on intuitions are not truly
challenging for state-of-the-art models. This situation is reminiscent of the
Turing test, which was believed to be a golden standard for proving machine
intelligence but is now considered insufficient for benchmarking the latest
language models [43]. The heuristically defined OOD tasks are akin to this
early benchmark, where existing models largely pass yet more challenging tests
are necessary to push the boundaries of what is considered state-of-the-art.
Choosing such tasks reflects mainly interpolation capabilities of ML models
but does not truly represent out-of-distribution generalizability. We advocate
for more attention and discussion on what constitutes a genuinely demanding
OOD task. This will pave the way for creating more convincing test scenarios,
with which we can identify, understand, and bridge the gaps towards truly
generalizable ML models.

##  IV Methods

###  IV.1 Datasets and models

We use the snapshots of the JARVIS, Materials Project, and OMQD databases used
in our previous study [38], available on Zenodo at
https://zenodo.org/records/8200972. The snapshots correspond the JARVIS
2022.12.12 version, the Materials Project 2021.11.10 version, and the OQMD
v1.6 version (published in November 2023), and they have been preprocessed to
remove materials with formation energies larger than 5 eV/atom.

Performing hyperparameter tuning for each of over 700 OOD tasks for every
model would be computationally impractical. In addition, hyperparameter tuning
techniques are developed mainly to optimize the in-distribution rather than
OOD performance. Our strategy is to perform the hyperparameter tuning on an
in-distribution validation set randomly sampled from the whole dataset; then
we apply the same optimal hyperparameters to every OOD task. Here, the
criterion for the optimal hyperparameters is based on not only model
performance but also training cost, to make the training amenable to extensive
benchmarking. Our objective is to demonstrate the overall generalization
performance on various OOD tasks rather than ensure optimal performance on
each OOD task for every model. We expect future work to develop techniques to
optimize OOD performance by focusing on a small set of representative tasks
selected based on this study.

For tree-based models, we use a descriptor set that consists of 145
compositional features [44] extracted from chemical formula and 128 structural
features [45] extracted from crystal structures. For the RF model [46], we use
100 estimators, 30 % of the features for the best splitting, and default
settings for other hyperparameters. For the XGB model [23], we use 4 parallel
boosting trees; for each boosting tree, we use 1000 estimators, a learning
rate of 0.25, an L1 (L2) regularization strength of 0.01 (0.1), and the
histogram tree grow method; we set the subsample ratio of training instances
to 0.85, the subsample ratio of columns to 0.3 when constructing each tree,
and the subsample ratio of columns to 0.5 for each level.

The deep learning models (ALIGNN, GMP, and LLM-Prop) are trained as follows.
ALIGNN models are trained with 2 ALIGNN layers, 2 GCN layers, a batch size of
16, 25 epochs, and layer normalization, while keeping other hyperparameters
the same as in the original implementation [26]. LLM-Prop models are trained
with a batch size of 64, a max token length of 1500, and a drop rate of 0.5,
while keeping the tokenizer, preprocessing strategy, and other hyperparameters
as in the original implementation [27]. For GMP models, we use the GMP-
featurizer package [47] to calculate the Gaussian Multipole features by using
40 radial probes, max 4th-order angular probes, and a Gaussian width of 1.5.
The Kresse-Joubert projector augmented-wave pseudopotentials are used to
calculated electron density [48]. The model architecture consists of a
compressor neural network (hidden layer dimension [512,256,128,64,32]), a
Set2Set neural network with 5 hidden layers to convert sets of varying sizes
to a fixed size tensor [49], and a predictor neural network (hidden layer
dimension [256,128,64,32,16]). GMP models are trained with 200 epochs, a batch
size of 128, and a ReduceLROnPlateau scheduler to tune the learning rate on a
validation set (10 % of the initial training set).

###  IV.2 SHAP analysis for leave-one-element-out tasks

The following procedure is proposed to identify the sources of biases in the
models trained for the leave-one-element-out tasks. For each task, we use the
model trained without element $X$ to make predictions on the whole dataset.
Then, we train a second model on the prediction errors of the first model on
the whole dataset. By doing so, the second model learns to correct the bias in
the first model which is the one used in the leave-one-element-out task.
Finally, we calculate the SHAP feature contributions of the second trained
model for all data with element $X$, interpreted as the correction for the
feature contribution of the first model for the out-of-distribution materials.
Here the first and second models are chosen to be XGB because it uses two
types of interpretable features: the compositional features are chemical
attributes computed from chemical formula [44], whereas the structural
features are characteristics of the local atomic environment calculated from
crystal structures [45]. We apply the above procedure to calculate the SHAP
contributions for compositional and structural features of the second XGB
model and aggregate the contributions into two categories.

###  IV.3 UMAP dimension reduction

The 2-component UMAP representation [39] are built using a local neighborhood
size of 50 and a minimum distance of 0.1. For deep learning models, the
embeddings are used as inputs to construct the UMAP plots. For tree-based
models, we sequentially drop highly correlated features with a Pearson
correlation threshold of 0.7 and then construct the UMAP plots based on the
final feature set. We use the kernel-density estimation method implemented in
SciPy [50] to calculate the probability density of training data for every
test data point. A probability density threshold of 0.001 is used to separate
between in-domain and out-of-domain test data.

## Data availability

We use the same snapshots as in our prior work [38] for the JARVIS, MP, and
OQMD datasets, which can be found on Zenodo at
https://zenodo.org/records/8200972.

## Code availability

The code used for ML training and analysis is available on GitHub at (to be
inserted upon acceptance).

## References

  * Butler  _et al._ [2018] K. T. Butler, D. W. Davies, H. Cartwright, O. Isayev, and A. Walsh, Machine learning for molecular and materials science, Nature 559, 547 (2018).
  * Carleo  _et al._ [2019] G. Carleo, I. Cirac, K. Cranmer, L. Daudet, M. Schuld, N. Tishby, L. Vogt-Maranto, and L. Zdeborová, Machine learning and the physical sciences, Reviews of Modern Physics 91, 045002 (2019).
  * Krenn  _et al._ [2022] M. Krenn, R. Pollice, S. Y. Guo, M. Aldeghi, A. Cervera-Lierta, P. Friederich, G. dos Passos Gomes, F. Häse, A. Jinich, A. Nigam, _et al._ , On scientific understanding with artificial intelligence, Nature Reviews Physics , 1 (2022).
  * Huang  _et al._ [2023] B. Huang, G. F. von Rudorff, and O. A. von Lilienfeld, The central role of density functional theory in the ai age, Science 381, 170 (2023).
  * Liu  _et al._ [2021] J. Liu, Z. Shen, Y. He, X. Zhang, R. Xu, H. Yu, and P. Cui, Towards out-of-distribution generalization: A survey, arXiv preprint arXiv:2108.13624 (2021).
  * Merchant  _et al._ [2023] A. Merchant, S. Batzner, S. S. Schoenholz, M. Aykol, G. Cheon, and E. D. Cubuk, Scaling deep learning for materials discovery, Nature , 1 (2023).
  * Yang  _et al._ [2024] H. Yang, C. Hu, Y. Zhou, X. Liu, Y. Shi, J. Li, G. Li, Z. Chen, S. Chen, C. Zeni, _et al._ , Mattersim: A deep learning atomistic model across elements, temperatures and pressures, arXiv preprint arXiv:2405.04967 (2024).
  * Batatia  _et al._ [2023] I. Batatia, P. Benner, Y. Chiang, A. M. Elena, D. P. Kovács, J. Riebesell, X. R. Advincula, M. Asta, W. J. Baldwin, N. Bernstein, _et al._ , A foundation model for atomistic materials chemistry, arXiv preprint arXiv:2401.00096 (2023).
  * Schmidt  _et al._ [2022] J. Schmidt, N. Hoffmann, H.-C. Wang, P. Borlido, P. J. Carriço, T. F. Cerqueira, S. Botti, and M. A. Marques, Large-scale machine-learning-assisted exploration of the whole materials space, arXiv preprint arXiv:2210.00579 (2022).
  * Deng  _et al._ [2023] B. Deng, P. Zhong, K. Jun, J. Riebesell, K. Han, C. J. Bartel, and G. Ceder, Chgnet as a pretrained universal neural network potential for charge-informed atomistic modelling, Nature Machine Intelligence 5, 1031 (2023).
  * Chen and Ong [2022] C. Chen and S. P. Ong, A universal graph deep learning interatomic potential for the periodic table, Nature Computational Science 2, 718 (2022).
  * Li  _et al._ [2024] K. Li, K. Choudhary, B. DeCost, M. Greenwood, and J. Hattrick-Simpers, Efficient first principles based modeling via machine learning: from simple representations to high entropy materials, Journal of Materials Chemistry A (2024).
  * Chen  _et al._ [2023] W. Chen, A. Hilhorst, G. Bokas, S. Gorsse, P. J. Jacques, and G. Hautier, A map of single-phase high-entropy alloys, Nature Communications 14, 2856 (2023).
  * Bokas  _et al._ [2021] G. B. Bokas, W. Chen, A. Hilhorst, P. J. Jacques, S. Gorsse, and G. Hautier, Unveiling the thermodynamic driving forces for high entropy alloys formation through big data ab initio analysis, Scripta Materialia 202, 114000 (2021).
  * Chen  _et al._ [2018] H.-L. Chen, H. Mao, and Q. Chen, Database development and calphad calculations for high entropy alloys: Challenges, strategies, and tips, Materials Chemistry and Physics 210, 279 (2018).
  * Frey  _et al._ [2023] N. C. Frey, R. Soklaski, S. Axelrod, S. Samsi, R. Gomez-Bombarelli, C. W. Coley, and V. Gadepally, Neural scaling of deep chemical models, Nature Machine Intelligence 5, 1297 (2023).
  * Kaplan  _et al._ [2020] J. Kaplan, S. McCandlish, T. Henighan, T. B. Brown, B. Chess, R. Child, S. Gray, A. Radford, J. Wu, and D. Amodei, Scaling laws for neural language models, arXiv preprint arXiv:2001.08361 (2020).
  * Choudhary  _et al._ [2020] K. Choudhary, K. F. Garrity, A. C. Reid, B. DeCost, A. J. Biacchi, A. R. Hight Walker, Z. Trautt, J. Hattrick-Simpers, A. G. Kusne, A. Centrone, _et al._ , The joint automated repository for various integrated simulations (jarvis) for data-driven materials design, npj computational materials 6, 173 (2020).
  * Wines  _et al._ [2023] D. Wines, R. Gurunathan, K. F. Garrity, B. DeCost, A. J. Biacchi, F. Tavazza, and K. Choudhary, Recent progress in the jarvis infrastructure for next-generation data-driven materials design, Applied Physics Reviews 10, 10.1063/5.0159299 (2023).
  * Jain  _et al._ [2013] A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, _et al._ , Commentary: The materials project: A materials genome approach to accelerating materials innovation, APL materials 1 (2013).
  * Saal  _et al._ [2013] J. E. Saal, S. Kirklin, M. Aykol, B. Meredig, and C. Wolverton, Materials design and discovery with high-throughput density functional theory: The open quantum materials database (OQMD), Jom 65, 1501 (2013).
  * Breiman [2001] L. Breiman, Random forests, Machine learning 45, 5 (2001).
  * Chen and Guestrin [2016] T. Chen and C. Guestrin, XGBoost, in  _Proc. 22nd ACM SIGKDD Int. Conf. Knowl. Discov. Data Min._ (ACM, New York, NY, USA, 2016) pp. 785–794.
  * Ward  _et al._ [2018] L. Ward, A. Dunn, A. Faghaninia, N. E. Zimmermann, S. Bajaj, Q. Wang, J. Montoya, J. Chen, K. Bystrom, M. Dylla, _et al._ , Matminer: An open source toolkit for materials data mining, Computational Materials Science 152, 60 (2018).
  * Lei and Medford [2022] X. Lei and A. J. Medford, A universal framework for featurization of atomistic systems, The Journal of Physical Chemistry Letters 13, 7911 (2022).
  * Choudhary and DeCost [2021] K. Choudhary and B. DeCost, Atomistic line graph neural network for improved materials property predictions, npj Computational Materials 7, 185 (2021).
  * Rubungo  _et al._ [2023] A. N. Rubungo, C. Arnold, B. P. Rand, and A. B. Dieng, Llm-prop: Predicting physical and electronic properties of crystalline solids from their text descriptions, arXiv preprint arXiv:2310.14029 (2023).
  * Riebesell  _et al._ [2022] J. Riebesell, R. Goodall, and S. G. Baird, Pymatviz: visualization toolkit for materials informatics (2022), https://github.com/janosh/pymatviz.
  * Deng  _et al._ [2024] B. Deng, Y. Choi, P. Zhong, J. Riebesell, S. Anand, Z. Li, K. Jun, K. A. Persson, and G. Ceder, Overcoming systematic softening in universal machine learning interatomic potentials by fine-tuning, arXiv preprint arXiv:2405.07105 (2024).
  * Choudhary and Sumpter [2023] K. Choudhary and B. G. Sumpter, Can a deep-learning model make fast predictions of vacancy formation in diverse materials?, AIP Advances 13 (2023).
  * Lundberg  _et al._ [2020] S. M. Lundberg, G. Erion, H. Chen, A. DeGrave, J. M. Prutkin, B. Nair, R. Katz, J. Himmelfarb, N. Bansal, and S.-I. Lee, From local explanations to global understanding with explainable ai for trees, Nature Machine Intelligence 2, 2522 (2020).
  * Behler and Parrinello [2007] J. Behler and M. Parrinello, Generalized neural-network representation of high-dimensional potential-energy surfaces, Physical review letters 98, 146401 (2007).
  * Li  _et al._ [2023a] K. Li, B. DeCost, K. Choudhary, M. Greenwood, and J. Hattrick-Simpers, A critical examination of robustness and generalizability of machine learning prediction of materials properties, npj Computational Materials 9, 55 (2023a).
  * Zhang  _et al._ [2023] H. Zhang, W. W. Chen, J. M. Rondinelli, and W. Chen, Et-al: Entropy-targeted active learning for bias mitigation in materials data, Applied Physics Reviews 10 (2023).
  * Schrier  _et al._ [2023] J. Schrier, A. J. Norquist, T. Buonassisi, and J. Brgoch, In pursuit of the exceptional: Research directions for machine learning in chemical and materials science, Journal of the American Chemical Society 145, 21699 (2023).
  * Goodall  _et al._ [2022] R. E. Goodall, A. S. Parackal, F. A. Faber, R. Armiento, and A. A. Lee, Rapid discovery of stable materials by coordinate-free coarse graining, Science Advances 8, eabn4117 (2022).
  * Ganose and Jain [2019] A. M. Ganose and A. Jain, Robocrystallographer: automated crystal structure text descriptions and analysis, MRS Communications 9, 874 (2019).
  * Li  _et al._ [2023b] K. Li, D. Persaud, K. Choudhary, B. DeCost, M. Greenwood, and J. Hattrick-Simpers, Exploiting redundancy in large materials datasets for efficient machine learning with less data, Nature Communications 14, 7283 (2023b).
  * McInnes  _et al._ [2018] L. McInnes, J. Healy, N. Saul, and L. Großberger, UMAP: Uniform Manifold Approximation and Projection, J. Open Source Softw. 3, 861 (2018).
  * Riebesell  _et al._ [2023] J. Riebesell, R. E. Goodall, A. Jain, P. Benner, K. A. Persson, and A. A. Lee, Matbench discovery–an evaluation framework for machine learning crystal stability prediction, arXiv preprint arXiv:2308.14920 (2023).
  * Faber  _et al._ [2017] F. A. Faber, L. Hutchison, B. Huang, J. Gilmer, S. S. Schoenholz, G. E. Dahl, O. Vinyals, S. Kearnes, P. F. Riley, and O. A. Von Lilienfeld, Prediction errors of molecular machine learning models lower than hybrid dft error, Journal of chemical theory and computation 13, 5255 (2017).
  * Wei  _et al._ [2022] J. Wei, Y. Tay, R. Bommasani, C. Raffel, B. Zoph, S. Borgeaud, D. Yogatama, M. Bosma, D. Zhou, D. Metzler, _et al._ , Emergent abilities of large language models, arXiv preprint arXiv:2206.07682 (2022).
  * Biever [2023] C. Biever, Chatgpt broke the turing test-the race is on for new ways to assess ai, Nature 619, 686 (2023).
  * Ward  _et al._ [2016] L. Ward, A. Agrawal, A. Choudhary, and C. Wolverton, A general-purpose machine learning framework for predicting properties of inorganic materials, npj Comput. Mater. 2, 1 (2016), 1606.09551 .
  * Ward  _et al._ [2017] L. Ward, R. Liu, A. Krishna, V. I. Hegde, A. Agrawal, A. Choudhary, and C. Wolverton, Including crystal structure attributes in machine learning models of formation energies via voronoi tessellations, Physical Review B 96, 024104 (2017).
  * Pedregosa  _et al._ [2011] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and É. Duchesnay, Scikit-learn: Machine Learning in Python, J. Mach. Learn. Res. 12, 2825 (2011).
  * Lei and Montoya [2023] X. Lei and J. Montoya, Gmp-featurizer: A parallelized python package for efficiently computing the gaussian multipole features of atomic systems, Journal of Open Source Software 8, 5476 (2023).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Physical review b 59, 1758 (1999).
  * Vinyals  _et al._ [2015] O. Vinyals, S. Bengio, and M. Kudlur, Order matters: Sequence to sequence for sets, arXiv preprint arXiv:1511.06391 (2015).
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, _et al._ , Scipy 1.0: fundamental algorithms for scientific computing in python, Nature methods 17, 261 (2020).

## Conflicts of interest

There are no conflicts of interest to declare.

Certain commercial products or company names are identified here to describe
our study adequately. Such identification is not intended to imply
recommendation or endorsement by the National Institute of Standards and
Technology, nor is it intended to imply that the products or names identified
are necessarily the best available for the purpose.

## Acknowledgments

This research was undertaken thanks in part to funding provided to the
University of Toronto’s Acceleration Consortium from the Canada First Research
Excellence Fund (Grant number: CFREF-2022-00042). Computational resources were
provided by the Calcul Quebec, Westgrid, and Compute Ontario consortia in the
Digital Research Alliance of Canada, and the Acceleration Consortium at the
University of Toronto. We thank Rhys Goodall for constructive comments and
discussion.

## Author contributions

K.L. and J.H.-S. conceived the idea. K.L. designed the project, performed the
investigation, and wrote the manuscript. A.N.R., X.L., and D.P. assisted in
the investigation. J.H.-S. supervised the project. All authors contributed to
the discussion and revision.