# Empirical Investigation of Neural Symbolic Reasoning Strategies

Yoichi Aoki,1 Keito Kudo,1 Tatsuki Kuribayashi,1,2 Ana Brassard,3,1
Masashi Yoshikawa,1 Keisuke Sakaguchi,1,3 Kentaro Inui1,3
1Tohoku University, 2Langsmith, Inc., 3RIKEN
<EMAIL_ADDRESS>
<EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS>
###### Abstract

Neural reasoning accuracy improves when generating intermediate reasoning
steps. However, the source of this improvement is yet unclear. Here, we
investigate and factorize the benefit of generating intermediate steps for
symbolic reasoning. Specifically, we decompose the reasoning strategy w.r.t.
step granularity and chaining strategy. With a purely symbolic numerical
reasoning dataset (e.g., A=1, B=3, C=A+3, C?), we found that the choice of
reasoning strategies significantly affects the performance, with the gap
becoming even larger as the extrapolation length becomes longer. Surprisingly,
we also found that certain configurations lead to nearly perfect performance,
even in the case of length extrapolation. Our results indicate the importance
of exploring effective strategies for neural reasoning models. 111Code
available at: https://github.com/ao1neko/reasoning-strategy

##  1 Introduction

Artificial intelligence researchers have been attempting neural-symbolic
integration for a long time d’Avila Garcez and Lamb (2020); Hamilton et al.
(2022). Neural models tend to perform better when generating intermediate
reasoning steps in addition to the answer. This phenomenon was seen across
various reasoning tasks, such as math word problems Wei et al. (2022); Cobbe
et al. (2021); Kojima et al. (2022); Recchia (2021); Lewkowycz et al. (2022),
commonsense reasoning Wei et al. (2022); Wang et al. (2022), and symbolic
reasoning Wei et al. (2022); Kojima et al. (2022). However, it is yet unclear
which factors in the intermediate step generation bring the benefit. Previous
studies often used different strategies for step generation in an ad-hoc
manner.

Figure 1: In a controlled setting, we found that output and chaining strategy
choice significantly impact performance when conducting multi-step reasoning.

(a) _All-at-once_ : output the entire reasoning chain and answer in a single
call. _Step-by-step_ : iteratively build the output with a single calculation
step per call. _Token-by-token_ : iteratively output only one _token_ per
call.

(b) The graph nodes represent variables and edges their dependencies.
_Shortest path_ : a minimal chain starting from the first necessary equation.
_Exhaustive_ : greedily solve all equations until the target is reached.
_Backward_ : start from the target’s equation, backtrack over dependencies
until a known value is reached, then solve each equation in order.

Figure 2: Overview of (a) output and (b) chaining strategies given the INPUT:
D=A+2, A=1, B=A+1, C=3+B, C?

To investigate this, we break down the neural reasoning process into two
strategies: _output strategy_ and _chaining strategy_ (Figure 1). The output
strategy (§2.1) determines the granularity of intermediate reasoning step
generation (all at once vs. step-by-step vs. token-by-token). Some studies
trained the models to generate reasoning steps and a conclusion derived from
them at once Nye et al. (2021); Lewkowycz et al. (2022); Wei et al. (2022);
Kojima et al. (2022); Wang et al. (2022); Recchia (2021), some generated a
single reasoning step given the input and iterated this process until
achieving a conclusion Sanyal et al. (2022); Picco et al. (2021); Tafjord et
al. (2021), and others iteratively generated sub-goals as well as reasoning
steps Liang et al. (2021); Shwartz et al. (2020).

In turn, the chaining strategy (§2.2) defines the reasoning path direction
(shortest path vs. exhaustive path vs. backward path). For example, some
studies used a backward chaining process Picco et al. (2021); Rocktäschel and
Riedel (2017); Cingillioglu and Russo (2019), while others adopted exhaustive
searches Tafjord et al. (2021); Liang et al. (2021); Yang et al. (2022).

To compare the strategies, we prepared a test bed of numerical reasoning
problems in a simplified language (Figure 1). This format allows for more
controlled testing while serving as a necessary condition—should a model fail
to solve it, it cannot be expected to adequately generalize to more complex
math word problems.

We found that both strategies substantially affect the symbolic reasoning
performance of neural seq2seq learners. Overall, iterative generation
outperformed all-at-once outputting, and roughly granular reasoning steps
(i.e., shortest-path chaining) lagged behind finely granular steps (i.e.,
exhaustive and backward chaining). Surprisingly, some settings had near-
perfect performance even in generalization tests which extrapolate over
greater reasoning depths and unseen numbers during training.

##  2 Experimental settings

#### Problem definition.

We evaluated the models’ ability to iteratively perform arithmetic operations
over given symbols. Given a series of equations, the task is to answer the
value of a target variable (Figure 1). Each question also has a certain
reasoning depth—the number of _necessary_ equations to reach the answer. For
example, the depth of the question A=1, B=2+A, C=3+B, D=2, C? is $3$ (A=1,
B=2+A, C=3+B).

Each equation defines either an assignment (e.g., A=1) or a modular addition
and an assignment (e.g., B=3+1). The addition is mod 100. The question
contexts also contain distractors that are not necessary to calculate the
answer (e.g., D=A+2 in Figure 1). A value assigned to a particular variable is
typically referred to in different equations (e.g., A=1, B=A+1). Numbers,
variables, and the ordering of equations are randomly assigned.

#### Motivation for using artificial data

There are mainly three advantages to this dataset. First, the symbolic format
allows easier control of reasoning depth for generalization tests.
Specifically, we trained a model using instances with shallow (1-5) depths and
evaluated them with instances with shallow/deep (1-12) depths. On the other
hand, math word problems are harder to control for reasoning depth (e.g., it
is not easy to come up with various instances which have a reasoning depth of
10). Second, we wanted to avoid the ”spurious bias” that natural (math word)
texts implicitly bring into the model Gururangan et al. (2018); Gupta et al.
(2021); Al-Negheimish et al. (2021); Sugawara et al. (2018); Jia and Liang
(2017); McCoy et al. (2019). Third, we assume that our setting is the
necessary condition for solving math word problems. It is unreasonable to
expect that a model that can’t solve this pure numerical reasoning task can
solve more complex tasks.

In total, we prepared 5K instances for training and 2.4K for testing.

###  2.1 Output strategies

We compared three configurations: all-at-once, step-by-step, and token-by-
token (Figure 2(a)).

All-at-once:  The model outputs the entire reasoning chain and the final
answer in a single call (i.e., chain-of-thought style) Wei et al. (2022);
Cobbe et al. (2021); Yavuz et al. (2022); Shwartz et al. (2020) . In this
setting, the more reasoning steps, the longer the sequence the decoder must
generate at once.

Step-by-step:  The model outputs a single reasoning step per call. Each
generated step is concatenated to the past input, and the model again
generates the next step (i.e., proofwriter style) Liang et al. (2021); Sanyal
et al. (2022); Picco et al. (2021); Tafjord et al. (2021); Shwartz et al.
(2020) . This process is iterated until the model outputs the answer or until
a set maximum number of iterations is reached ($100$). Token-by-token:  This
is the same as step-be-step chaining, but the decoder outputs only a single
_token_ per call. We set the maximum number of steps to $500$.

Comparing all-at-once and the others reveals the effect of changing the
sequence length that the decoder outputs in a single call. In addition,
comparing step-by-step and token-by-token quantifies the advantage of breaking
a problem into meaningful units.

###  2.2 Chaining strategies

Particular variables sometimes depend on another variable; the key to reaching
the correct answer is determining the order in which the equations are
referred to. Regarding existing studies, we compared three chaining
strategies: _shortest-path_ , _exhaustive_ , and _backward_ chaining (Figure
2(b)).

Shortest-path chaining:  The model straight-forwardly solves the equations
starting from the first solvable one (i.e., involving a known value) and
ending with the target Wei et al. (2022); Cobbe et al. (2021); Yavuz et al.
(2022); Shwartz et al. (2020). Here, the reasoning behind determining the
shortest path is not outputted by the model.

Exhaustive chaining:  The model greedily solves all given equations until the
target value is reached Tafjord et al. (2021); Liang et al. (2021); Yang et
al. (2022). Specifically, the model calculates the left-most solvable equation
in each step. Note that this strategy typically derives a long reasoning
chain; from an engineering perspective, this strategy is inefficient.

Backward chaining:  The model starts from the equation for the target variable
and backtracks over the dependent equations until it reaches a known value
Picco et al. (2021); Rocktäschel and Riedel (2017); Cingillioglu and Russo
(2019). Then, it solves each equation in order by inserting known or
calculated values until the target one is reached.

No chaining:  As a baseline, we also examined the setting where the model was
trained to directly output the answer.

Figure 3: Distributions of the total reasoning chain length (num. characters).
The all-at-once and step-by-step generate those at depth 12.

(a) Output strategy

(b) Chaining strategy

Figure 4: Accuracy changes of the models against reasoning depth. The gray
range represents the training data domain (1-5 depth). Figure (4(a)) shows the
performance degradation with the increase of reasoning steps when using the
all-at-once strategy. Figure (4(b)) shows that the combination of step-by-step
output and backward/exhaustive chaining leads to successful generalization.

##  3 Results

#### Models:

We used the pre-trained T5-base, T5-large
222https://huggingface.co/docs/transformers/model_doc/T5 Raffel et al. (2020),
and BART-base 333https://huggingface.co/docs/transformers/model_doc/bart Lewis
et al. (2020). Results of BART-base are in Appendix C.

Note that their pre-defined tokenizers have all the numbers from 0 to 9, and
the numerical values in our dataset are divided into digits (e.g., “12” should
be “@@1 @@2”) in advance, following Kim et al. (2021).

Training:  The models were first pre-trained using a 10K simple dataset for 30
epochs, then trained with the 5K training set (1K training instances for each
reasoning depth.) for 2000 epochs. The experiment setting details are in
Appendix A. In addition, we prepared 0.2K test instances for each reasoning
depth. This pre-training is intended to teach the models primitive operations
(i.e., assignment, reference, and addition). The pre-training dataset contains
two types of single-depth instances: assign-refer type (e.g., A=1,A?) and
operate-assign-refer type (e.g., A=1+3, A?). All the results in the paper are
averages of the results on three different seeds.

###  3.1 Output strategies

We compared the output strategies while fixing the chaining strategy to the
shortest path. Figure 4(a) shows the accuracy per reasoning depth. Note that
the accuracy score here denotes whether the answer (e.g., C=6) is correct. We
observed the following: (i) generating intermediate reasoning steps enhance
the performance, and (ii) among the output strategies, step-by-step works the
best, and all-at-once works the worst. The format of the dataset in this study
is simple. Therefore, this result indicates the low symbolic reasoning ability
of neural models and the necessity of the choice of an appropriate reasoning
strategy.

We hypothesized that the source of all-at-once’s inferiority was that the
decoder overfitted to output a similar length of reasoning steps as those in
the (shallower) training data. In fact, the models generated relatively
shorter reasoning steps in the out-of-domain (e.g., depth of 12) setting when
using the all-at-once strategy (Figure 3); this supports our hypothesis.

The advantage of step-by-step over token-by-token suggests the advantage of
breaking the problem into meaningful units (reasoning step) and modeling each
step in a single call of the encoder-decoder.

Depth | Shortest | Backward | Exhaustive
---|---|---|---
6 | 99.3/99.3 | 100/ 100 | 99.7/99.7
8 | 95.5/95.7 | 100/ 100 | 99.8/99.8
12 | 76.7/77.7 | 99.5/99.5 | 98.2/98.3
Table 1: Accuracy of the T5-base model with the step-by-step output strategy
at each depth (chain/answer). Question: A=1, B=2+A, B?
---
Error types |  Gold |  Prediction
Copying error |  B=2+A, |  B=6+A,
|  B=2+1, |  B=6+1,
|  B=3 |  B=7
Hasty assignment |  B=2+A, |  (skip step)
|  B=2+1, |  B=2+2,
|  B=3 |  B=4
Table 2: Illustrative examples of the errors under the step-by-step, shortest-
path chaining settings. (skip step) denotes that the reasoning steps is
accidentally skipped.

###  3.2 Chaining strategies

Figure 4(b) and Table 1 show the results on each depth with a fixed step-by-
step output strategy. Note that the accuracy of the chain (left side of the
scores) was measured based on not an exact match but mathematically. For
example, even if the order of generated equations is different, it is correct.
The results of a fixed token-by-token output strategy are in Appendix B.

While the performance dropped in the shortest-path setting as the reasoning
depth increased, with either the exhaustive or backward chaining, models
successfully solved the task even when extrapolating to depths 6-12. The
models correctly generated the intermediate steps (nearly perfect) as well as
the final answer in the exhaustive and backward chaining settings (Table 1).
Note that these strategies were ineffective with all-at-once outputting.

Gontier et al. (2020) compared chaining strategies and concluded that models
that didn’t generate reasoning steps had better generalization performance
than models that did when the reasoning chains were long. However, our results
suggest that the choice of the appropriate output strategy improves the
reasoning ability of the model.

We considered that the source of shortest-path inferiority was the rough
granularity of the given reasoning steps. The models don’t know the shortest
path before outputting the reasoning steps. Therefore, both the exhaustive and
shortest path chaining approaches must search for variables other than those
on the shortest path. As shown in Figure 2(b), the exhaustive chaining
approach is taught this process explicitly. On the other hand, the shortest-
path chaining approach must be learned that by training data that don’t
include this process. We thought this difference affected the accuracy and
concluded that the accuracy is higher when the granularity of given
intermediate steps is finer, even though they are long.

Therefore, we concluded that the accuracy is higher when the granularity of
intermediate steps is finer, even though they are long.

###  3.3 Error analysis

We also analyzed the errors of the depth-12 instances under the shortest-path
strategy. 444In total, 32 instances were analyzed. That is the total number of
incorrect answers on one seed. We observed two types of errors: (i) copying
errors and (ii) hasty assignment. Table 2 shows an illustrative example of
each error type and the percentage of these errors. The most frequent one
(53%) was a simple copying error, where the model failed to accurately copy an
original equation into the reasoning chain. This erroneous copying ability is
consistent with Xu et al. (2020) and supports the advantage of introducing a
copy mechanism to the model Ontanon et al. (2022). Second, a hasty assignment
is the model skipping the step to copy the equation from context and instead
assigned it a random value. Note that these errors were almost addressed in
the other strategies; this could stem from the difficulty of the implicit
calculation of the shortest path.

###  3.4 Models’ scalability

To investigate the scalability, we compared T5-large with T5-base. Figure 5
shows the result. T5-large had a similar trend but slightly lower accuracy on
all-at-once and step-by-step compared to T5-base. The reason may be that
T5-large needs more data for updating the weights of the entire model. On the
other hand, the accuracy of T5-large is higher than T5-base on token-by-token.
It’s because the data size of token-by-token is as token lengths of output
sequence times as the data size of all-at-once, as shown in Figure 2(a). This
result indicates that the parameter size of the model needs to be larger to
output token-by-token.

Figure 5: Accuracy changes of the T5-base and T5-large against reasoning
depth. The gray range presents the training data domain (1-5 depth). This
figure shows that the accuracy of T5-large with token-by-token is higher.

##  4 Conclusions

We investigated and factorized the reasoning strategy in symbolic numerical
reasoning with neural seq2seq models. We found that the combination of step-
by-step output and finely granular reasoning leads to successfully performing
symbolic reasoning. Our results support the potential of neural models for
symbolic reasoning.

## Limitations

We found that even simple symbolic reasoning requires the appropriate
selection of reasoning strategy. It is unclear whether our findings generalize
to more complex symbolic reasoning and/or problems written in natural
language. If our findings do not generalize in these different settings, we
must address the gap in future work. For example, we start with one of the
simplest tasks and find out when models fail as we add complexity to tasks one
by one.

From the engineering perspective, the iterative strategies are limited to the
input length of the model. For example, in our experiments, when adopting the
setting where reasoning depths are greater than 13, the input length of step-
by-step and token-by-token became longer than the input length limit of T5
(i.e., 512 tokens).

In addition, gigantic language models (e.g., GPT-3) have recently been used.
Including these models in our study is one of our future works.

## Acknowledgements

We thank four anonymous reviewers who provided valuable feedback. We would
like to also appreciate the member of Tohoku NLP Group for their cooperation
in conducting this research.

This work was supported by JSPS KAKENHI Grant Numbers JP22H00524, 21K21343 and
JST CREST Grant Number JPMJCR20D2, Japan.

## References

  * Al-Negheimish et al. (2021) Hadeel Al-Negheimish, Pranava Madhyastha, and Alessandra Russo. 2021.  Numerical reasoning in machine reading comprehension tasks: are we there yet? In _Proceedings of the 2021 Conference on Empirical Methods in Natural Language Processing, EMNLP 2021, Virtual Event / Punta Cana, Dominican Republic, 7-11 November, 2021_ , pages 9643–9649. Association for Computational Linguistics.
  * Cingillioglu and Russo (2019) Nuri Cingillioglu and Alessandra Russo. 2019.  Deeplogic: Towards end-to-end differentiable logical reasoning.  In _Proceedings of the AAAI 2019 Spring Symposium on Combining Machine Learning with Knowledge Engineering (AAAI-MAKE 2019) Stanford University, Palo Alto, California, USA, March 25-27, 2019., Stanford University, Palo Alto, California, USA, March 25-27, 2019_ , volume 2350 of _CEUR Workshop Proceedings_. CEUR-WS.org.
  * Cobbe et al. (2021) Karl Cobbe, Vineet Kosaraju, Mohammad Bavarian, Jacob Hilton, Reiichiro Nakano, Christopher Hesse, and John Schulman. 2021.  Training verifiers to solve math word problems.  _CoRR_ , abs/2110.14168.
  * d’Avila Garcez and Lamb (2020) Artur d’Avila Garcez and Luís C. Lamb. 2020.  Neurosymbolic AI: the 3rd wave.  _CoRR_ , abs/2012.05876.
  * Gontier et al. (2020) Nicolas Gontier, Koustuv Sinha, Siva Reddy, and Christopher Pal. 2020.  Measuring systematic generalization in neural proof generation with transformers.  In _Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems 2020, NeurIPS 2020, December 6-12, 2020, virtual_.
  * Gupta et al. (2021) Ashim Gupta, Giorgi Kvernadze, and Vivek Srikumar. 2021.  BERT & family eat word salad: Experiments with text understanding.  In _Thirty-Fifth AAAI Conference on Artificial Intelligence, AAAI 2021, Thirty-Third Conference on Innovative Applications of Artificial Intelligence, IAAI 2021, The Eleventh Symposium on Educational Advances in Artificial Intelligence, EAAI 2021, Virtual Event, February 2-9, 2021_ , pages 12946–12954. AAAI Press.
  * Gururangan et al. (2018) Suchin Gururangan, Swabha Swayamdipta, Omer Levy, Roy Schwartz, Samuel R. Bowman, and Noah A. Smith. 2018.  Annotation artifacts in natural language inference data.  In _Proceedings of the 2018 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, NAACL-HLT, New Orleans, Louisiana, USA, June 1-6, 2018, Volume 2 (Short Papers)_ , pages 107–112. Association for Computational Linguistics.
  * Hamilton et al. (2022) Kyle Hamilton, Aparna Nayak, Bojan Bozic, and Luca Longo. 2022.  Is neuro-symbolic AI meeting its promise in natural language processing? A structured review.  _CoRR_ , abs/2202.12205.
  * Jia and Liang (2017) Robin Jia and Percy Liang. 2017.  Adversarial examples for evaluating reading comprehension systems.  In _Proceedings of the 2017 Conference on Empirical Methods in Natural Language Processing, EMNLP 2017, Copenhagen, Denmark, September 9-11, 2017_ , pages 2021–2031. Association for Computational Linguistics.
  * Kim et al. (2021) Jeonghwan Kim, Giwon Hong, Kyung-min Kim, Junmo Kang, and Sung-Hyon Myaeng. 2021\.  Have you seen that number? investigating extrapolation in question answering models.  In _Proceedings of the 2021 Conference on Empirical Methods in Natural Language Processing, EMNLP 2021, Virtual Event / Punta Cana, Dominican Republic, 7-11 November, 2021_ , pages 7031–7037. Association for Computational Linguistics.
  * Kojima et al. (2022) Takeshi Kojima, Shixiang Shane Gu, Machel Reid, Yutaka Matsuo, and Yusuke Iwasawa. 2022.  Large language models are zero-shot reasoners.  _CoRR_ , abs/2205.11916.
  * Lewis et al. (2020) Mike Lewis, Yinhan Liu, Naman Goyal, Marjan Ghazvininejad, Abdelrahman Mohamed, Omer Levy, Veselin Stoyanov, and Luke Zettlemoyer. 2020.  BART: denoising sequence-to-sequence pre-training for natural language generation, translation, and comprehension.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics, ACL 2020, Online, July 5-10, 2020_ , pages 7871–7880. Association for Computational Linguistics.
  * Lewkowycz et al. (2022) Aitor Lewkowycz, Anders Andreassen, David Dohan, Ethan Dyer, Henryk Michalewski, Vinay V. Ramasesh, Ambrose Slone, Cem Anil, Imanol Schlag, Theo Gutman-Solo, Yuhuai Wu, Behnam Neyshabur, Guy Gur-Ari, and Vedant Misra. 2022\.  Solving quantitative reasoning problems with language models.  _CoRR_ , abs/2206.14858.
  * Liang et al. (2021) Zhengzhong Liang, Steven Bethard, and Mihai Surdeanu. 2021.  Explainable multi-hop verbal reasoning through internal monologue.  In _Proceedings of the 2021 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, NAACL-HLT 2021, Online, June 6-11, 2021_ , pages 1225–1250. Association for Computational Linguistics.
  * McCoy et al. (2019) Tom McCoy, Ellie Pavlick, and Tal Linzen. 2019.  Right for the wrong reasons: Diagnosing syntactic heuristics in natural language inference.  In _Proceedings of the 57th Conference of the Association for Computational Linguistics, ACL 2019, Florence, Italy, July 28- August 2, 2019, Volume 1: Long Papers_ , pages 3428–3448. Association for Computational Linguistics.
  * Nye et al. (2021) Maxwell I. Nye, Anders Johan Andreassen, Guy Gur-Ari, Henryk Michalewski, Jacob Austin, David Bieber, David Dohan, Aitor Lewkowycz, Maarten Bosma, David Luan, Charles Sutton, and Augustus Odena. 2021.  Show your work: Scratchpads for intermediate computation with language models.  _CoRR_ , abs/2112.00114.
  * Ontanon et al. (2022) Santiago Ontanon, Joshua Ainslie, Zachary Fisher, and Vaclav Cvicek. 2022.  Making transformers solve compositional tasks.  In _Proceedings of the 60th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 3591–3607, Dublin, Ireland. Association for Computational Linguistics.
  * Picco et al. (2021) Gabriele Picco, Thanh Lam Hoang, Marco Luca Sbodio, and Vanessa López. 2021\.  Neural unification for logic reasoning over natural language.  In _Findings of the Association for Computational Linguistics: EMNLP 2021, Virtual Event / Punta Cana, Dominican Republic, 16-20 November, 2021_ , pages 3939–3950. Association for Computational Linguistics.
  * Raffel et al. (2020) Colin Raffel, Noam Shazeer, Adam Roberts, Katherine Lee, Sharan Narang, Michael Matena, Yanqi Zhou, Wei Li, and Peter J. Liu. 2020.  Exploring the limits of transfer learning with a unified text-to-text transformer.  _J. Mach. Learn. Res._ , 21:140:1–140:67.
  * Recchia (2021) Gabriel Recchia. 2021.  Teaching autoregressive language models complex tasks by demonstration.  _CoRR_ , abs/2109.02102.
  * Rocktäschel and Riedel (2017) Tim Rocktäschel and Sebastian Riedel. 2017.  End-to-end differentiable proving.  In _Advances in Neural Information Processing Systems 30: Annual Conference on Neural Information Processing Systems 2017, December 4-9, 2017, Long Beach, CA, USA_ , pages 3788–3800.
  * Sanyal et al. (2022) Soumya Sanyal, Harman Singh, and Xiang Ren. 2022.  Fairr: Faithful and robust deductive reasoning over natural language.  In _Proceedings of the 60th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers), ACL 2022, Dublin, Ireland, May 22-27, 2022_ , pages 1075–1093. Association for Computational Linguistics.
  * Shwartz et al. (2020) Vered Shwartz, Peter West, Ronan Le Bras, Chandra Bhagavatula, and Yejin Choi. 2020\.  Unsupervised commonsense question answering with self-talk.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing, EMNLP 2020, Online, November 16-20, 2020_ , pages 4615–4629. Association for Computational Linguistics.
  * Sugawara et al. (2018) Saku Sugawara, Kentaro Inui, Satoshi Sekine, and Akiko Aizawa. 2018.  What makes reading comprehension questions easier? In _Proceedings of the 2018 Conference on Empirical Methods in Natural Language Processing, Brussels, Belgium, October 31 - November 4, 2018_ , pages 4208–4219. Association for Computational Linguistics.
  * Tafjord et al. (2021) Oyvind Tafjord, Bhavana Dalvi, and Peter Clark. 2021.  Proofwriter: Generating implications, proofs, and abductive statements over natural language.  In _Findings of the Association for Computational Linguistics: ACL/IJCNLP 2021, Online Event, August 1-6, 2021_ , volume ACL/IJCNLP 2021 of _Findings of ACL_ , pages 3621–3634. Association for Computational Linguistics.
  * Wang et al. (2022) Xuezhi Wang, Jason Wei, Dale Schuurmans, Quoc V. Le, Ed H. Chi, and Denny Zhou. 2022\.  Self-consistency improves chain of thought reasoning in language models.  _CoRR_ , abs/2203.11171.
  * Wei et al. (2022) Jason Wei, Xuezhi Wang, Dale Schuurmans, Maarten Bosma, Ed H. Chi, Quoc Le, and Denny Zhou. 2022.  Chain of thought prompting elicits reasoning in large language models.  _CoRR_ , abs/2201.11903.
  * Xu et al. (2020) Song Xu, Haoran Li, Peng Yuan, Youzheng Wu, Xiaodong He, and Bowen Zhou. 2020.  Self-attention guided copy mechanism for abstractive summarization.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics, ACL 2020, Online, July 5-10, 2020_ , pages 1355–1362. Association for Computational Linguistics.
  * Yang et al. (2022) Kaiyu Yang, Jia Deng, and Danqi Chen. 2022.  Generating natural language proofs with verifier-guided search.  _CoRR_ , abs/2205.12443.
  * Yavuz et al. (2022) Semih Yavuz, Kazuma Hashimoto, Yingbo Zhou, Nitish Shirish Keskar, and Caiming Xiong. 2022.  Modeling multi-hop question answering as single sequence prediction.  In _Proceedings of the 60th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers), ACL 2022, Dublin, Ireland, May 22-27, 2022_ , pages 974–990. Association for Computational Linguistics.

Figure 6: Accuracy changes of token-by-token per reasoning depth. The gray
range presents the training data domain (depths 1-5). Figure 7: Accuracy
changes of the T5-base and BART-base models per reasoning depth. The gray
range presents the training data domain (depths 1-5). T5 seems to outperform
BART.

##  Appendix A Details on Experimental Settings

We first examined the learning rate from $10^{-3}$, $10^{-4}$, and $10^{-5}$;
among them, we used the largest rate at which the loss converged. After
training models, we used the model with the lowest validation loss among the
per-epoch checkpoints during the training reported. We used four NVIDIA V100
GPUs for NVLink 16GiB HBM2.

##  Appendix B Results of Token-by-token

Figure 6 shows the results on each depth with a fixed token-by-token output
strategy. Like step-by-step, the performance drops in the shortest-path
setting as the reasoning depth increases. In addition, the exhaustive or
backward successfully solves the task even when extrapolating to depths 6-12.

##  Appendix C Different Architectures

We also tested BART-base Lewis et al. (2020) as a baseline to investigate the
effectiveness of the NLP-task-oriented objectives used in the T5-style pre-
training. Figure 7 shows this result. In this particular setting, T5 was
superior to BART. This suggests that the NLP-task-oriented objectives benefit
symbolic reasoning.

Question: A=1, C=5+B, B=2+A, D=3+A, C?
---
Chain error types |  Gold |  Prediction
Ignoring the incorrect step |  A=1, B=2+A, B=2+1, B=3, C=5+B, C=5+3, C=8 |  A=1, B=2+D, B=2+A, B=2+1, B=3, C=5+B, C=5+3, C=8
Correct assignment |  A=1, B=2+A, B=2+1, B=3, C=5+B, C=5+3, C=8 |  A=1, B=2+D, B=2+1, B=3, C=5+B, C=5+3, C=8
Non affecting error |  A=1, B=2+A, B=2+1, B=3, C=5+B, C=5+3, C=8 |  A=1, B=2+A, B=2+1, B=3, D=3+A, D=3+2, D=5, C=5+B, C=5+3, C=8
Table 3: These instances are examples of chain errors. Note that the final
answers are correct.

##  Appendix D Other errors

We analyzed the cases where the answer is correct and the chain is wrong.
Table 3 shows examples of chain errors. Ignoring the incorrect step is an
example of the model outputting the correct reasoning step after outputting an
incorrect one. Correct assignment is an example in which the assignment
accidentally makes the model output the correct step. Finally, Non-affecting
error is an example in which a variable not on the shortest path is wrongly
assigned a value.


[1,2]RaghulParthipan [3]Hannah M.Christensen [2,4]J. ScottHosking [1]Damon
J.Wischik

1]Department of Computer Science, University of Cambridge, Cambridge, UK
2]British Antarctic Survey, Cambridge, UK 3]Department of Physics, University
of Oxford, Oxford, UK 4]The Alan Turing Institute, London, UK

Raghul Parthipan<EMAIL_ADDRESS>
# Using Probabilistic Machine Learning to Better Model Temporal Patterns in
Parameterizations: a case study with the Lorenz 96 model.

###### Abstract

The modelling of small-scale processes is a major source of error in climate
models, hindering the accuracy of low-cost models which must approximate such
processes through parameterization. Red noise is essential to many operational
parameterization schemes, helping model temporal correlations. We show how to
build on the successes of red noise by combining the known benefits of
stochasticity with machine learning. This is done using a physically-informed
recurrent neural network within a probabilistic framework. Our model is
competitive and often superior to both a bespoke baseline and an existing
probabilistic machine learning approach (GAN) when applied to the Lorenz 96
atmospheric simulation. This is due to its superior ability to model temporal
patterns compared to standard first-order autoregressive schemes. It also
generalises to unseen scenarios. We evaluate across a number of metrics from
the literature, and also discuss the benefits of using the probabilistic
metric of hold-out likelihood.

A major source of inaccuracies in climate models is due to ‘unresolved’
processes. These occur at scales smaller than the resolution of the climate
model (‘sub-grid’ scales) but still have key effects on the overall climate.
In fact, most of the inter-model spread in how much global surface
temperatures increase after $\textrm{CO}{{}_{2}}$ concentrations double is due
to the representation of clouds (Schneider et al., 2017; Zelinka et al.,
2020). The typical approach to deal with the problem of unresolved processes
has been to model the effects of these unresolved processes as a function of
the resolved ones. This is known as ‘parameterization’.

Red noise is a key feature in many parameterization schemes (Christensen et
al., 2015; Johnson et al., 2019; Molteni et al., 2011; Palmer et al., 2009;
Skamarock et al., 2019; Walters et al., 2019). Hasselmann (1976) developed the
theory underpinning this, showing that the coupling of processes with
different time-scales leads to red noise signals in the longer time-scale
process, analogous to what takes place in Brownian motion. In
parameterization, the resolved processes typically have far longer time-scales
than the unresolved ones, motivating the inclusion of red noise in these
schemes. The usefulness of tracking the history of the system for
parameterization follows from the importance of red noise.

We use probabilistic machine learning to propose a data-driven successor to
red noise in parameterization. The theory explaining the prevalence of red
noise in the climate was developed in an idealized model, within the framework
of physics-based differential equations. Such approaches may be intractable
outside of the idealized case. Recently, there has been much work looking at
uncovering relationships from data instead (Arcomano et al., 2022; Beucler et
al., 2020, 2021; Bolton and Zanna, 2019; Brenowitz and Bretherton, 2018, 2019;
Chattopadhyay et al., 2020a; Gagne et al., 2020; Gentine et al., 2018;
Krasnopolsky et al., 2013; O’Gorman and Dwyer, 2018; Rasp et al., 2018;
Vlachas et al., 2022; Yuval and O’Gorman, 2020; Yuval et al., 2021). We
develop a Recurrent Neural Network in a probabilistic framework, combining the
benefits of stochasticity with the ability to learn temporal patterns in more
flexible ways than permitted by red noise. We use the Lorenz 96 model (Lorenz,
1996), henceforth the L96, as a proof-of-concept. Our model is competitive
with, and often outperforms, a bespoke baseline. It also generalises to unseen
scenarios.

###  0.1 Numerical Models

Numerical models represent the state of the Earth system at time $t$ by a
state vector $\mathbf{X}_{t}$. The components of $\mathbf{X}_{t}$ may include,
for example, the mean temperature and humidity at time $t$ in every cell of a
grid that covers the Earth. The goal is to parameterize the effects of
unresolved processes in such a way that the model can reproduce the evolution
of this finite-dimensional representation of the state of the real Earth
system. A simple way to model the evolution would be

$\displaystyle X_{k,t+1}$ $\displaystyle=X_{k,t}+f_{k}(\mathbf{X}_{t})$ (1)

where $X_{k,t}$ is the value that the state variable $X$ takes at the spatial
coordinate $k$ and time point $t$ and $X_{k,t}\in\mathbb{R}^{d}$,
$\mathbf{X}_{t}\in\mathbb{R}^{dK}$, and $f$ is a function for the updating
process.

###  0.2 Stochastic Schemes

Stochasticity can be included using the following form

$\displaystyle X_{k,t+1}$
$\displaystyle=X_{k,t}+f_{k}(\mathbf{X}_{t},h_{k,t+1})$ (2) $\displaystyle
h_{k,t+1}$ $\displaystyle=\beta(h_{k,t},z_{k,t+1})$ (3) $\displaystyle
z_{k,t}$ $\displaystyle\sim\mathcal{N}(0,I)$ (4)

where invented hidden variables (discussed below) are denoted
$h_{k,t}\in\mathbb{R}^{H}$ and are tracked through time, $\beta$ is a function
for updating these, and $z_{k,t}\in\mathbb{R}^{Z}$ is a source of
stochasticity.

Introducing stochasticity through (2) improved on models of the form (1).
Results included better ensemble forecasts (Buizza et al., 1999; Leutbecher et
al., 2017; Palmer, 2012), and improvements to model mean state (Berner et al.,
2012) and climate variability (Christensen et al., 2017). The motivation for
using stochasticity comes from the understanding that the effects of the
unresolved (sub-grid) processes cannot be effectively predicted as a
deterministic function of the resolved ones due to a lack of scale separation
between them. Stochasticity allows us to capture our uncertainty about those
aspects of the unresolved processes which may affect the resolved outcomes.

###  0.3 Hidden Variables and Red Noise

Hidden variables — $h_{k,t}$ in (2) and (3) — are defined here as being
variables separate to the observed state, $X_{k,t}$, which if tracked help
better model $X_{k,t}$. Using them is key to numerical weather and climate
models using stochastic parameterizations, allowing temporal correlations to
be better modelled. Red noise results when the hidden variables evolve with a
first-order autoregressive (AR1) process such as

$h_{k,t+1}=\phi h_{k,t}+\sigma(1-\phi^{2})^{1/2}z_{k,t+1}$ (5)

One example is the stochastically perturbed parameterization tendencies (SPPT)
scheme (Buizza et al., 1999; Palmer et al., 2009) which is widely used in
forecasting models (Leutbecher et al., 2017; Molteni et al., 2011; Palmer et
al., 2009; Sanchez et al., 2016; Skamarock et al., 2019; Stockdale et al.,
2011). Here, the AR1 process results in far better weather and climate
forecasting skill than a simple white noise model, with good modelling of
regime behaviour requiring correlated noise (Christensen et al., 2015; Dawson
and Palmer, 2015).

There is no intrinsic reason why an AR1 process is the best way to deal with
these correlations. It is simply a modelling choice.

###  0.4 Machine Learning for Parameterization

Learning the parameters of a climate model, either the simple form (1) or the
general form (2)–(4), requires deciding on a parametric form for the functions
$f$ and $\beta$. For full-scale climate models, it is difficult to find
appropriate functions. The machine learning (ML) approach is to _learn_ these
from data. Various researchers have proposed ML methods for learning the
deterministic model (1) (Beucler et al., 2020, 2021; Bolton and Zanna, 2019;
Brenowitz and Bretherton, 2018, 2019; Gentine et al., 2018; Krasnopolsky et
al., 2013; O’Gorman and Dwyer, 2018; Rasp et al., 2018; Yuval and O’Gorman,
2020; Yuval et al., 2021).

Amongst ML-trained stochastic models (Gagne et al., 2020; Guillaumin and
Zanna, 2021), various ones with red noise were proposed by Gagne et al.
(2020), using Generative Adversarial Networks (GANs) (Goodfellow et al.,
2014). Full details of the architecture are in Gagne et al. (2020) and we will
refer to one of their best-performing models (which they call X-sml-r) as the
_GAN_. A wider range of generative models can be trained using such an
advesarial approach as opposed to maximum likelihood (the standard way to
train models in ML, discussed in Section 2) but such methods are notoriously
unstable due to the nature of the minimax loss in training. Although they used
ML to learn $f$ in (2), it was not used to learn $\beta$, which they modelled
with an AR1 process.

Recurrent Neural Networks (RNNs) are a popular ML tool for modelling
temporally correlated data, eliminating the need for update functions to be
manually specified. RNNs have had great success in the ML literature in a
variety of sequence modelling tasks, including text generation (Graves, 2013;
Sutskever et al., 2011), machine translation (Sutskever et al., 2014) and
music generation (Eck and Schmidhuber, 2002; Mogren, 2016). The state-of-the-
art RNNs are gated ones, principally the long short-term memory (LSTM)
networks (Hochreiter and Schmidhuber, 1997) and gated recurrent unit (GRU)
networks (Cho et al., 2014), which provide major improvements to the standard
issue of unstable gradients.

Current parameterization work using ML to model temporal correlations has used
deterministic approaches, including deterministic RNNs and echo state networks
(Arcomano et al., 2022; Chattopadhyay et al., 2020a, b; Vlachas et al., 2018,
2020, 2022) and found notable success. This suggests such ML approaches are
effective ways to model temporal dynamics in physical systems. However, these
studies do not incorporate stochasticity.

Other off-the-shelf models are not obviously suited for the parameterization
task. Transformers and attention-based models (Vaswani et al., 2017) perform
well for sequences but require all the previous data to be tracked for each
simulation step, providing a computational burden which increases simulation
cost. Random Forests (RFs) have been used for parameterization (O’Gorman and
Dwyer, 2018; Yuval and O’Gorman, 2020) and shown to be stable at run-time (due
to predicting averages form the training set) but it is not obvious how they
would learn and track hidden variables.

###  0.5 Overview

The L96 set-up and baselines are presented in Section 1. Our model is detailed
in Section 2, with the experiments and results in Section 3. This is followed
by a discussion on the use of ‘likelihood’ in evaluating probabilistic climate
models (Section 4), with our conclusions in Section 2.

##  1 Parameterization in the Lorenz 96

We introduce the L96 model here and then present two L96 parameterization
models from the literature which help clarify the above discussion and serve
as our baselines.

###  1.1 Lorenz 96 Model

We use the two-tier L96 model, a toy model for atmospheric circulation that is
extensively used for stochastic parameterization studies (Arnold et al., 2013;
Crommelin and Vanden-Eijnden, 2008; Gagne et al., 2020; Kwasniok, 2012; Rasp,
2020). We use the configuration described in Gagne et al. (2020). It consists
of two scales of variables: a large, low-frequency variable, $X_{k}$, coupled
to small, high-frequency variables, $Y_{j,k}$. These are dimensionless
quantities, evolving as follows:

$\frac{dX_{k}}{dt}=\underbrace{-X_{k-1}(X_{k-2}-X_{k+1})}_{\textrm{advection}}\underbrace{-X_{k}}_{\textrm{diffusion}}\
\underbrace{+F}_{\textrm{forcing}}-\underbrace{\frac{hc}{b}\sum_{j=J(k-1)+1}^{kJ}Y_{j}}_{\textrm{coupling}}\mkern
70.0muk=1,...,K$
$\frac{dY_{j,k}}{dt}=\underbrace{-cbY_{j+1,k}(Y_{j+2,k}-Y_{j-1,k})}_{\textrm{advection}}\>\underbrace{-\,cY_{j,k}}_{\textrm{diffusion}}\>\underbrace{-\frac{hc}{b}X_{k}}_{\textrm{coupling}}\mkern
150.0muj=1,...,J$

where in our experiments, the number of $X_{k}$ variables is $K=8$, and the
number of $Y_{j,k}$ variables per $X_{k}$ is $J=32$. The value of the
constants are set to $h=1,\>b=10$ and $\>c=10$. These indicate that the fast
variable evolves ten times faster than the slow variable and has one-tenth the
amplitude.

The L96 model is good for parameterization work as we can consider the $X_{k}$
to be coarse processes resolved in both low-resolution and high-resolution
simulators, whilst the $Y_{j,k}$ can be considered as those that can only be
resolved in high-resolution, computationally expensive simulators. In this
study the coupled set of L96 equations are treated as the ‘truth’, and the aim
is to learn a good model for the evolution of $X$ alone using this ‘truth’
data. The effects of $Y_{j,k}$ must therefore be parameterized. Success would
be if the modelled evolution of $X$ matched that from the truth.

The L96 is also useful as it contains separate persistent dynamical regimes,
which change with different values of $F$ (Christensen et al., 2015; Lorenz,
2006). The dominant regime exhibits a wave-two pattern around the ring of $X$
variables, whilst the rarer regime exhibits a wave-one type pattern. In the
atmosphere, regimes include persistent circulation patterns like the
Pacific/North American (PNA) pattern and the North Atlantic Oscillation (NAO).

It is easy to use models with no memory which do well on standard loss metrics
but fail to capture this interesting regime behaviour. This is seen in the
L96, where including temporally correlated (red) noise improves the statistics
describing regime persistence and frequency of occurrence (Christensen et al.,
2015) as well as improving weather forecasting skill (Arnold et al., 2013;
Palmer, 2012). Temporally correlated noise also improves regime behaviour
(frequency of occurrence and persistence) in operational models (Dawson and
Palmer, 2015), as well as forecasting skill (Palmer et al., 2009).

###  1.2 Parameterization Models in the Literature

The below models are stochastic. They use AR1 processes to model temporal
correlations, but this need not be the case, as we show in Section 2.

####  1.2.1 Stochastic Non-ML Model with Non-ML Hidden Variables

Arnold et al. (2013) propose a model which we refer to as the _polynomial_
model (in light of the polynomial in (6)) which is

$\displaystyle X_{k,t+1}$
$\displaystyle=X_{k,t}+\omega_{k}(\mathbf{X}_{t})-\Delta
t\,\Bigl{(}aX_{k,t}^{3}+bX_{k,t}^{2}+cX_{k,t}+d+h_{k,t+1}\Bigr{)}$ (6)

where $h_{k,t}$ evolves as in (5) with $h_{k,1}=\sigma z_{k,1}$, and where

$\omega_{k}(\mathbf{X}_{t})=\lambda_{k}\Big{(}\mathbf{X}_{t}+\lambda(\mathbf{X}_{t}/2)\Big{)}$
(7)

and

$\lambda_{k}(\mathbf{X}_{t})=\Delta
t\,\Big{(}-X_{k-1,t}(X_{k-2,t}-X_{k+1,t})-X_{k,t}+F\Big{)}$ (8)

where (7) is the implementation of a second order Runge-Kutta method and
$[\lambda(\mathbf{a})]_{k}=\lambda_{k}(\mathbf{a})$. $h_{k,t}$ only depends on
$h_{k,t-1}$.

####  1.2.2 Stochastic ML Model with Non-ML Hidden Variables

The GAN from Gagne et al. (2020) replaces (6) by

$\displaystyle X_{k,t+1}$
$\displaystyle=X_{k,t}+\omega_{k}(\mathbf{X}_{t})-\Delta t\
U(X_{k,t},h_{k,t+1})$ (9)

where the function $U$ is implemented by a neural network (NN), $\omega$ is
defined in (7)–(8), $h_{k,t}$ evolves as in (5) with $\sigma=1$, and
$h_{k,1}=z_{k,1}$.

##  2 Our Proposed RNN Model

###  2.1 Stochastic ML Model with ML Hidden Variables

Our model, henceforth denoted the RNN, follows the general form in (2)–(4) but
splits the hidden variable into two parts: $h_{k,t}=(r_{k,t},l_{k,t})$. The
model is

$\displaystyle X_{k,t+1}$
$\displaystyle=X_{k,t}+\omega_{k}(\mathbf{X}_{t})-\Delta
t\Big{(}g_{\theta}(X_{k,t})+r_{k,t+1}\Big{)}$ (10) $\displaystyle r_{k,t+1}$
$\displaystyle=b_{\theta}(l_{k,t+1})+\sigma z_{k,t+1}$ (11) $\displaystyle
l_{k,t+1}$ $\displaystyle=s_{\theta}(l_{k,t},r_{k,t})$ (12)

where $g$, $b$ and $s$ are NNs with weights $\theta$, $z_{k,t}$ is exogenous
noise as defined in (4), $\omega$ is specified in equation (7) above
(implements a second order Runge-Kutta method and expressing physical
quantities like advection), and $\theta$ and $\sigma$ are parameters to be
learnt. The dimensions of the NNs are:
$g_{\theta}:\mathbb{R}^{1}\rightarrow\mathbb{R}^{1}$,
$b_{\theta}:\mathbb{R}^{8}\rightarrow\mathbb{R}^{1}$ and
$s_{\theta}:\mathbb{R}^{9}\rightarrow\mathbb{R}^{8}$. It is structurally based
on a Recurrent Neural Network (RNN) as it has hidden variables, the same
update procedure is used each step, and $g$, $b$ and $s$ are NNs.

The key insight is that our model allows more flexibility than the standard
AR1 processes for expressing temporal correlations. To clarify this, equation
(10) is identical in form to (6), but the hidden state is evolved in a more
flexible manner than the AR1 process in (5).

Figure 1 shows the mechanism of generation and the NN architectures used to
learn the functions. The dotted lines denote the action of deterministic
functions, with the solid lines showing the sources of stochasticity. The main
architectural details are that $s$ in (12) is implemented using two GRU
layers, each composed of four units, and $b$ in (11) is represented with a
dense layer of size one. $g$ in (10) is implemented using three fully-
connected layers.

Here, as well as in the baselines, the parameterization models are ‘spatially
local’ meaning that the full $\mathbf{X}_{t}$ vector is not taken in as input
when modelling $X_{k,t+1}$ anywhere apart from in
$\omega_{k}(\mathbf{X}_{t})$. This is done to mimic parameterizations in
operational weather and climate models. For all forecast models, $\Delta
t=0.005$ model time units (MTU). One MTU is approximately five atmospheric
days when considering the time it takes for errors to double.

Figure 1: RNN graphical model showing how each $\mathbf{X}_{t}$ is generated.
$\mathbf{X}_{t+1}$ is a function of $\mathbf{X}_{t}$ and $\mathbf{r}_{t+1}$,
and $\mathbf{r}_{t}$ is a function of $\mathbf{l}_{t}$ and $\mathbf{z}_{t}$.
$\mathbf{z}_{t}$ is the source of stochasticity. $b_{\theta}$, $s_{\theta}$
and $g_{\theta}$ are NNs. $\mathbf{l}_{t}$ consists of a stack of $l_{t}^{1}$
and $l_{t}^{2}$, each $\in\mathbb{R}^{4}$.

###  2.2 Training Probability Models Using Likelihood

The RNN is probabilistic and trained using likelihood. The ‘likelihood of a
sequence’ is denoted $\mathrm{Pr}(x_{1},x_{2},...,x_{n})$ and can be
interpreted as the probability assigned to a given sequence of variables
$(x_{1},x_{2},...,x_{n})$ by a given model.

For our RNN, the log-likelihood (the natural logarithm of the likelihood) of
the sequence of $\mathbf{X}_{t}$ is

$\log\mathrm{Pr}(\mathbf{x}_{1},...,\mathbf{x}_{n};\mathbf{x}_{0})=\sum_{k=1}^{K}\Big{(}\log\mathrm{Pr}(r_{k,1})+\sum_{t=2}^{n}\log\mathrm{Pr}(r_{k,t}|l_{k,t})\Big{)}-Kn\log(\Delta
t)$

where $r$ and $l$ are deterministic functions of
$\mathbf{x}_{0},...,\mathbf{x}_{n}$ derived from equations (10)–(12). The full
derivation is in A.

A common way to train probabilistic models is to maximise the likelihood of
the training data. This involves finding a set of parameters ($\theta$ and
$\sigma$ here) which maximise the likelihood. Other enhancements to training
involve including regularisation penalties such as dropout. The explicit form
of the RNN’s likelihood makes training relatively easy.

###  2.3 RNN Training

Truth data for training was created by running the full L96 model five
separate times with the following values of $F$: $(19,20,20.5,21,21.5)$. We
gave our model data from different forcing scenarios so it could learn how
perturbations in the forcing may affect the evolution of $\mathbf{X}_{t}$. The
GAN and polynomial were trained only on $F=20$ as was done by their authors.
The truth data was created by solving the two-level L96 equations using a
fourth-order Runge-Kutta timestepping scheme and a fixed time step $\Delta
t=0.001$ MTU, with the output saved at every 0.005 MTU. A training set of
length 2,500 MTU was assembled from the truth data, consisting of three equal
500 MTU components with $F=(19,20.5,21)$ and one 1,000 MTU component with
$F=20$, with a $F=21.5$ set of length 500 MTU kept as a validation hold-out
set.

The RNN was trained using truncated back propagation through time on sequences
of length 700 time steps for 100 epochs with a batch size of 32 using Adam
(Kingma and Ba, 2014), with $\theta$ and $\sigma$ being the learnable
parameters. A variable learning rate was used, starting at 0.0001 for the
first 70 epochs and decayed to 0.00003 for the remaining 30. The model
parameters which gave the lowest loss on the validation set were saved.

##  3 Results

This section analyses performance across a range of time-scales. The first
results are for $F=20$. Assessing the model in the training realm allows us to
verify if it can replicate the L96 attractor. Later experiments are for $F\geq
28$. For $F=20$ and $F=28$, 50,000 MTU long simulations ($\approx 685$
‘atmospheric years’) were created for analysis.

###  3.1 Weather Evaluation

The models were evaluated in a weather forecast framework. 745 initial
conditions were randomly selected from the truth attractor and an ensemble of
40 forecasts each lasting 3.5 MTU were produced from each initial condition.
Figure 2 shows the spread and error terms for these experiments over time. The
error is defined as

$\text{error}(t)=\sqrt{\frac{1}{M}\sum_{m=1}^{M}\big{(}X_{m}^{O}(t)-\overline{X_{m}^{\text{sample}}(t)}\big{)}^{2}}$

where there are $M$ different initial conditions, $X_{m}^{O}(t)$ is the
observed state at time $t$ for the $m^{\text{th}}$ initial condition, and
$\overline{X_{m}^{\text{sample}}(t)}$ is the ensemble mean forecast at time
$t$, initialized at $t_{init}$, such that $t=t_{init}+\tau$.

The spread is defined as

$\text{spread}(t)=\sqrt{\frac{1}{M}\sum_{m=1}^{M}\frac{1}{N}\sum_{n=1}^{N}\big{(}X_{m,n}(t)-\overline{X_{m}^{\text{sample}}(t)}\big{)}^{2}}$

where $N$ is the number of ensemble members and $X_{m,n}(t)$ is the state of
the $n^{\text{th}}$ member at time $t$ for the $m^{\text{th}}$ initial
condition.

As noted by Leutbecher and Palmer (2008), for a perfectly reliable forecast
(defined as one where $X_{m,n}(t)$ and $X_{m}^{O}(t)$ are independent samples
from the same distribution), for large $N$ and $M$ the error should be equal
to the spread. A reduction in error indicates an ensemble forecast that better
tracks observations, and a spread/error ratio close to one indicates a
reliable forecast. Figure 2 shows that the polynomial and RNN have the
smallest errors, but the RNN has the better spread/error ratio after 0.75 MTU.
The GAN is underdispersive, with the greatest error but smallest spread/error
ratio, suggesting it is overconfident in its predictions.

Figure 2: Error and spread for weather experiments. We expect a smaller error
for a more accurate forecast model. The spread/error ratio would be close to
one in a ‘perfectly reliable’ forecast.

###  3.2 Climate Evaluation

The ability to simulate the climate of the L96 was evaluated using the 50,000
MTU simulations. The histograms in Figure 3a are approximations of each
model’s marginal distribution of $X_{k,t}$. Successful reproduction of the L96
climate would result in a model’s histogram matching the truth model’s.
Qualitatively, all the models are similar to the truth.

Figure 3: Histograms of $X_{k,t}$. (a) Comparisons for the full simulation
data. (b) Depiction which shows sampling error. The truth histogram is shown
in green and bold. The closer a model’s histogram to the truth, the better it
models the probability density of $X_{k,t}$.

We can quantitatively compare how well histograms match the truth using KL-
divergence. We compute this as follows: we discretize our modelled continuous
distribution of a variable by binning into a histogram, and do the same for
the truth model. We then evaluate the KL-divergence between $q(a)$, the
distribution described by the histogram of a variable $A$ in the true L96
model, and $p(a)$, the distribution described by the histogram of that
variable in our model

$\text{KL}(q(a)||p(a))=\sum_{a}q(a)\log\frac{q(a)}{p(a)}$

where the sum is over all the discrete values which $A$ takes. The smaller
this measure, the better the goodness-of-fit. The results (Table 1 column one)
confirm that the RNN best matches the truth model here.

Table 1: KL-divergence (goodness-of-fit) between 1) the truth and 2) the distributions shown in the noted Figures. The smaller the KL-divergence, the better the match between the true and modelled distributions. The best model in each case is shown in bold. Model | Figure 3a | Figure 4 | Figure 5a | Figure 5b | Figure 7a | Figure 7b
---|---|---|---|---|---|---
Polynomial | 0.004 | 3.9 | 0.03 | 0.040 | - | -
RNN | 0.001 | 1.8 | 0.02 | 0.004 | 0.008 | 0.006
GAN | 0.010 | 11.6 | 0.14 | 0.050 | 0.009 | 0.004

It it is also important to check whether the histograms are composed of enough
data such that they are good representations of each model’s actual
distribution of $X_{k,t}$. Figure 3b seeks to answer this by accounting for
sampling error. By splitting a particular model’s simulation run into five
equally sized sets and plotting these histograms on top of each other,
variability can be made apparent. Having done this, there is no noticeable
difference between Figure 3a and b, suggesting that the error from the
sampling process is small relative to the differences between each model.

The L96 model used here displays two distinct regimes with separate dynamics.
We use the approach from Christensen et al. (2015) to examine the regimes.
First, the time series are temporally smoothed with a running average over 0.4
MTU to help identify regimes (Stephenson et al., 2004; Straus et al., 2007).
The dimensionality of the system is then reduced using Principal Component
Analysis, as often done when studying atmospheric data (Selten and Branstator,
2004; Straus et al., 2007). For the truth time series, the components $PC1$
and $PC2$ are degenerate and are in phase quadrature, representing wavenumber
two oscillations. $PC3$ and $PC4$ are also degenerate and in phase quadrature,
representing wavenumber one oscillations. All model runs are projected onto
the truth model’s components. Given the degeneracies, the magnitude of the
principal component vectors, $||[PC1,PC2]||$ and $||[PC3,PC4]||$ are computed
and histograms of the system are plotted in this space.

Figure 4: Regime characteristics of the L96. 2D histograms show the magnitudes
of the projections of $X_{k,t}$ on to the principle components from the truth
series.

The presence of two regimes is apparent in Figure 4a where there is a large
maximum corresponding to the major regime, A, but also a smaller peak
corresponding to the minor regime, B. The polynomial fails to capture the two
regimes, whereas the ML models succeed. The RNN puts an appropriate amount of
density in Regime B, whilst the GAN puts too much. Comparison is made easier
by decomposing the 2D Figure 4 into two, 1D density plots (Figure 5). We can
use the KL-divergence to measure how well our models match the truth across
all three histograms. Table 1 shows the RNN best matches the truth and so
captures the regime characteristics best. There are still deficiencies though
as it does not put enough density in the right-hand side tail of Figure 5a.

Figure 5: Density plots for each regime for $F=20$. The truth is shown in
green and bold. (a) Magnitude of [$PC1$,$PC2$]. (b) Magnitude of
[$PC3$,$PC4$].

###  3.3 Generalisation Experiments

We set the L96 forcing to $F=28,32,35$ and $40$, and examine how the models
can capture changes due to varying external forcings. These $F$ values are
notably different to those in the training data so allow us to test
generalisation. For example, the change from $F=21.5$ (validation set) to
$F=28$ results in the following changes to the regime structure: the centroid
location of the rarer regime shifts to higher values of $||[PC1,PC2]||$ and
$||[PC3,PC4]||$ (Figure 6), and the proportion of time spent in the rarer
regime increases from $38\%$ to $50\%$. The range of $X_{k,t}$ also increases
from $[-12,19]$ to $[-24,29]$. The differences in wave behaviour between
regimes means the above changes result in different system dynamics.

Figure 6: Difference between the $F=21.5$ and $F=28$ densities in principal
component space. The principal components are those from the $F=20$ truth
series.

In all the below experiments the polynomial model exploded so is omitted. This
is merely due to the specification of a third-order polynomial. It
significantly deviates from its target values when $X_{k,t}$ values notably
different from those in training (such as $X_{k,t}\geq 19$) are taken as
inputs.

The to simulate the $F=28$ climate was explored in a similar manner to the
$F=20$ one. First, for the histograms of $X_{k,t}$ (analogous to Figure 3) the
RNN and GAN both had small KL divergences (0.0015 vs 0.0025). Next, the
principal component projections were examined (analogous to Figure 4). The
same components as determined from the $F=20$ truth data set are used. In this
space, the RNN has a smaller KL-divergence (0.6) than the GAN (0.7). Figures
for the above two plots are omitted due to it being hard to visually
distinguish the models. As before, further comparison can be made by examining
the two, 1D density plots in Figure 7 (analogous to Figure 5). Both models
perform well (KL divergence in Table 1), but neither put enough density in the
right-hand side tails.

Figure 7: Density plots for each regime for $F=28$. The truth is shown in
green and bold. (a) Magnitude of [$PC1$,$PC2$]. (b) Magnitude of
[$PC3$,$PC4$]. The right-hand side tails are not giving sufficient density by
the RNN nor the GAN.

The models were also evaluated in the weather forecasting framework for
$F=28,32,35$ and $40$. 750 initial conditions were randomly selected from the
truth attractors and an ensemble of 40 forecasts each lasting 2.5 MTU were
produced. Figure 8 shows the RNN and GAN have similar errors, but the RNN’s
spread is best matched with its error, whereas the GAN continues to be
underdispersive.

Figure 8: Error and spread for weather experiments for varying forcing values.
The RNN’s spread is best matched to its error unlike the GAN.

###  3.4 Computational Cost

We wish our models to have a lower cost than the full L96. The computational
cost is calculated by considering the number of floating-point operations per
simulated time step ($\Delta t=0.005$). The RNN (8,682) is notably cheaper
than the truth model (88,000) and the GAN (14,074), so meets the computational
cost objective of parameterization work. The polynomial (334) is the cheapest.

##  4 Evaluation and Diagnostics Using Likelihood

###  4.1 Hold-out Likelihood

We also evaluate using the likelihood (explained in Section 2.2) of hold-out
data. This is a standard approach in the ML literature. Evaluation is done on
hold-out sets to ensure that overly complex models which just overfit the
training data are not selected.

Here, hold-out sets for $F=20$ and $F=28$ were created by taking a 10,000 MTU
subset from the 50,000 MTU truth sets created in Section 3. The hold-out log-
likelihoods are shown in Table 2. The likelihood for the polynomial model and
its derivation is in A. We also present an approach to approximate the GAN
likelihood. This is despite its full form being intractable (involving
integrals which cannot be efficiently approximated using Monte Carlo sampling)
and so it not being typically used to evaluate GANs. This is also detailed in
A. The RNN has a lower likelihood on $F=28$ than the polynomial, despite the
polynomial exploding. This is discussed below.

###  4.2 What is the Use of Likelihood

####  4.2.1 For Evaluation

Other metrics only capture snapshots of model performance. Likelihood is a
composite measure which assesses a model’s full joint distribution. For
example, in the univariate case, the mean-squared-error and variance of a
model tell two separate things about its performance, with implicit
assumptions being made about variables being normally distributed (whenever
mean-squared-error is used). The likelihood captures both of these, and
without such assumptions.

Likelihood also captures information about more complex, joint distributions,
saving the need for custom metrics to be invented to assess specific features.
To illustrate this, consider we wish to assess a model’s temporal correlation.
The likelihood already contains this information. Although custom metrics
could be invented to assess this, the likelihood is an off-the-shelf metric
which is already available. We suggest it is wasteful not to use it.

Being a composite metric brings challenges though. Poor performance in certain
aspects may be overshadowed by good performance elsewhere. This can result in
cases where increased hold-out likelihoods do not correspond to better sample
quality, as noted in the ML literature (Goodfellow et al., 2014; Grover et
al., 2018; Theis et al., 2015; Zhao et al., 2020) and seen in our results: the
RNN has a worse $F=28$ hold-out likelihood than the polynomial, yet the
polynomial model explodes unlike the RNN. In this case, the phenomenon of
‘explosion’ is not significantly penalising the likelihood. Likelihood is
therefore a complement, not replacement, of other metrics which are important
to the end-user.

####  4.2.2 For Diagnostics

Likelihood’s composite nature makes it a helpful diagnostic tool. Just like
KL-divergence, the further away a model is — in any manner — from the data in
the hold-out set, the worse the likelihood will be. If a model has a poor
likelihood despite performing well on a range of standard metrics, this
suggests there are still deficiencies in the model which need investigating.
For example, with the RNN at $F=28$, the poor average likelihood is caused by
a tiny number of segments of the $X_{k,t}$ sequences being extremely poorly
modelled (and therefore having extremely poor likelihoods assigned). On
inspection, the issue is due to the choice of a fixed $\sigma$ in equation
(11) — this is appropriate for most of the time series, but for parts which
are difficult to model it is too small, preventing the model from expressing
sufficient uncertainty. This could be rectified by allowing $\sigma$ to vary,
and we suggest future work explore this.

Table 2: Log-likelihood on hold-out data. The RNN has the best likelihood in both cases. Model | $F=20$ | $F=28$
---|---|---
Polynomial | $4.98$ | 4.05
RNN | 6.53 | $2.69$
GAN | $\approx-37.69\pm 0.01$ | $\approx-147.02\pm 0.08$

[Discussion and Conclusion]

We present an approach to replace red noise with a more flexible stochastic
machine-learnt model for temporal correlations. Even though we used ML to
model $g_{\theta}(X_{k,t})$ (often called the ‘sub-grid tendency’) in (10),
the real benefit comes from using ML to model the hidden variables. For
example, on setting $g_{\theta}(X_{k,t})$ in our model equal to
$aX_{k,t}^{3}+bX_{k,t}^{2}+cX_{k,t}+d$ from the polynomial baseline, our model
performance hardly suffered (and the cost halved) for $F=20$. This points to a
physics-based ML approach where conventional parameterizations are augmented
with a stochastic term learnt using ML.

Using physical knowledge to structure ML models can help with learning. We
leveraged many elements from the polynomial baseline, particularly the the
physical features, and this gave us better results than when we modelled the
system without them. Although in theory a NN can _learn_ to create helpful
features (such as advection), giving useful features can make learning easier.
The NN does not need to learn the known, useful physical relationships and
instead can focus on learning other helpful ones.

There are more sophisticated ways to model the system. We noticed that in some
cases the RNN struggled even to overfit the training data, regardless of the
RNN complexity. This could be due to difficulties in learning the evolution of
the hidden variables. Creating architectures that permit better hidden
variables to be learnt (ones which model long-term correlations better) would
give better models. Despite our use of GRU cells, which along with the LSTM
were great achievements in sequence modelling, we still faced issues with the
modelling of long-term trends. Certain models (mainly those which did not
leverage physical structure) resulted in unstable simulations. Using
hierarchical models which learn to model trends at different temporal
resolutions (Chung et al., 2016; Liu et al., 2020) may provide improvement.

Learning from all the high-resolution data whilst still being computationally
efficient at simulation time is an interesting idea. Whilst we cannot keep all
the high-resolution dimensions in our schemes (as otherwise we might as well
run the high-resolution model outright — which is not possible as it is too
slow for the desired timescales and is a reason why this field of work
exists), only keeping the average (as in existing work which coarse-grains
data) means much data is thrown away. Our preliminary investigations used the
graphical model in Figure 9. We designed this so that in training,
$\mathbf{l}_{t}$, is learnt such that it contains useful information about the
high-resolution data, $\mathbf{Y}_{t}$, whereas during generation
$\mathbf{Y}_{t}$ is not required to be simulated. For the L96, this showed no
improvement over our RNN though.

Figure 9: Model which allows learning from high-resolution data without
requiring it at simulation time, where $N_{t}\sim\mathcal{N}(0,I)$.

Further work with GANs could result in better models. There may be aspects of
realism which we either may not be aware of or may not be able to quantify, so
are difficult to include explicitly in our probability models. The GAN
discriminator could learn what features constitute ‘realism’ and so the
generator may learn to create sequences which contain these. However, our work
using GANs and Wasserstein-GANs (Arjovsky et al., 2017; Gulrajani et al.,
2017) to train our RNN (instead of by likelihood) did not show benefit. This
may be due to the difficulty of training GANs on longer sequences. They have
been used to train RNNs on medical timeseries (Esteban et al., 2017) and to
model music (Mogren, 2016) so applying this approach to climatic sequences may
be something for future work to reconcile. The challenges of using GANs are
seen by how we have not been able to reproduce the results shown in Gagne et
al. (2020) despite the same set-up. This points to the instability of GAN
training as noted by them too. Mode collapse is a common issue encountered
when training GANs, where the generator fails to produce samples that explore
certain modes of the distribution, and could explain why Regime B in Figure 4
has more density than desired.

###  4.3 Conclusion

We have shown how to build on the benefits of red noise in parameterizations
by using a probabilistic ML approach based on an RNN. This now needs testing
on more complex systems — both to specifically improve on AR1 processes, and
to learn new, more flexible models of $\mathbf{X}_{t}$. An interesting area
for further work is to examine what information is tracked in the hidden
variables of the RNN, and how this relates to the timescales of the required
memory in the modelled system. The field is ripe for other probabilistic ML
tools to be used and we suspect that further customisation of these will lead
to many improved parameterization models.

Likelihood can often be fairly easily calculated, and where this be the case,
we propose that the community also evaluate the hold-out likelihood for any
devised probabilistic model (ML or otherwise). It is a useful debugging tool,
assessing the full joint distribution of a model. It would also provide a
consistent evaluation metric across the literature. Given the challenges
relating to sample quality, likelihood should complement not replace existing
metrics.

Finally, we have used demanding tests to show that ML models can generalise to
unseen scenarios. We cannot hope for ML models to generalise to all settings.
We do not expect this from our physics-based models either. But ML models
_can_ generalise outside of their training realm. And it is by using
challenging hold-out tests that we can assess their ability to do so, and if
they fail, begin the diagnosis of why.

All code used in this study (including the code required to create the data)
is publicly available at https://github.com/raghul-parthipan/l96_rnn.

##  Appendix A Likelihood Derivations

In all cases, the log-likelihood of the sequence of $\mathbf{X}_{t}$ is given
by

$\displaystyle\begin{split}\log\mathrm{Pr}(\mathbf{x}_{1},...,\mathbf{x}_{n};\mathbf{x}_{0})&=\log\mathrm{Pr}(\mathbf{x}_{1};\mathbf{x}_{0})+\log\mathrm{Pr}(\mathbf{x}_{2}|\mathbf{x}_{1};\mathbf{x}_{0})+\log\mathrm{Pr}(\mathbf{x}_{3}|\mathbf{x}_{2},\mathbf{x}_{1};\mathbf{x}_{0})\\\
{}&\qquad\qquad+...+\log\mathrm{Pr}(\mathbf{x}_{n}|\mathbf{x}_{n-1},...,\mathbf{x}_{1};\mathbf{x}_{0})\end{split}$
(13)

which follows from the laws of probability.

###  A.1 RNN

Now from (10), we can write

$\displaystyle(\mathbf{X}_{t}|\mathbf{X}_{t-1}{=}\mathbf{x}_{t-1},...,\mathbf{X}_{1}{=}\mathbf{x}_{1};\mathbf{X}_{0}{=}\mathbf{x}_{0})$
$\displaystyle=f(\mathbf{x}_{t-1})-\Delta
t(\mathbf{r}_{t}|\mathbf{X}_{t-1}{=}\mathbf{x}_{t-1},...,\mathbf{X}_{1}{=}\mathbf{x}_{1};\mathbf{X}_{0}{=}\mathbf{x}_{0})$
(14)

and given this relationship between $\mathbf{X}_{t}$ and $\mathbf{r}_{t}$, a
change of variables for the likelihood gives

$\log\mathrm{Pr}(\mathbf{x}_{t}|\mathbf{x}_{t-1},...,\mathbf{x}_{1};\mathbf{x}_{0})=\log\mathrm{Pr}(\mathbf{r}_{t}|\mathbf{r}_{t-1},...,\mathbf{r}_{1})-K\log(\Delta
t)$ (15)

(13) is therefore

$\displaystyle\begin{split}&=\log\mathrm{Pr}(\mathbf{r}_{1})+\log\mathrm{Pr}(\mathbf{r}_{2}|\mathbf{r}_{1})+\log\mathrm{Pr}(\mathbf{r}_{3}|\mathbf{r}_{2},\mathbf{r}_{1})\\\
{}&\qquad\qquad+...+\log\mathrm{Pr}(\mathbf{r}_{n}|\mathbf{r}_{n-1:1})-Kn\log(\Delta
t)\end{split}$ (16)
$\displaystyle=\log\mathrm{Pr}(\mathbf{r}_{1})+\sum_{t=2}^{n}\log\mathrm{Pr}(\mathbf{r}_{t}|\mathbf{l}_{t})-Kn\log(\Delta
t)$ (17)
$\displaystyle=\sum_{k=1}^{K}\Big{(}\log\mathrm{Pr}(r_{k,1})+\sum_{t=2}^{n}\log\mathrm{Pr}(r_{k,t}|l_{k,t})\Big{)}-Kn\log(\Delta
t)$ (18)

where (17)–(18) follow from the independencies in the graphical model (Figure
1).

###  A.2 Polynomial Model

For the polynomial model, the log-likelihood of the $\mathbf{X}_{t}$ sequence
is

$\displaystyle\log\mathrm{Pr}(\mathbf{x}_{1},...,\mathbf{x}_{n};\mathbf{x}_{0})$
$\displaystyle=\sum_{k=1}^{K}\Big{(}\log\mathrm{Pr}(h_{k,1})+\sum_{t=1}^{n-1}\log\mathrm{Pr}(h_{k,t+1}|h_{k,t})\Big{)}-Kn\log(\Delta
t)$ (19)

and the derivation follows below.

The graphical model is shown in Figure 10a. From (6), we can write

$\displaystyle(\mathbf{X}_{t}|\mathbf{X}_{t-1}{=}\mathbf{x}_{t-1},...,\mathbf{X}_{1}{=}\mathbf{x}_{1};\mathbf{X}_{0}{=}\mathbf{x}_{0})$
$\displaystyle=f(\mathbf{x}_{t-1})-\Delta
t(\mathbf{h}_{t}|\mathbf{X}_{t-1}{=}\mathbf{x}_{t-1},...,\mathbf{X}_{1}{=}\mathbf{x}_{1};\mathbf{X}_{0}{=}\mathbf{x}_{0})$
(20)

and given this relationship between $\mathbf{X}_{t}$ and $\mathbf{h}_{t}$, a
change of variables for the likelihood gives

$\log\mathrm{Pr}(\mathbf{x}_{t}|\mathbf{x}_{t-1},...,\mathbf{x}_{1};\mathbf{x}_{0})=\log\mathrm{Pr}(\mathbf{h}_{t}|\mathbf{h}_{t-1},...,\mathbf{h}_{1})-K\log(\Delta
t)$ (21)

(13) is therefore

$\displaystyle\begin{split}&=\log\mathrm{Pr}(\mathbf{h}_{1})+\log\mathrm{Pr}(\mathbf{h}_{2}|\mathbf{h}_{1})+\log\mathrm{Pr}(\mathbf{h}_{3}|\mathbf{h}_{2},\mathbf{h}_{1})\\\
{}&\qquad\qquad+...+\log\mathrm{Pr}(\mathbf{h}_{n}|\mathbf{h}_{n-1:1})-Kn\log(\Delta
t)\end{split}$ (22)
$\displaystyle=\log\mathrm{Pr}(\mathbf{h}_{1})+\sum_{t=1}^{n-1}\log\mathrm{Pr}(\mathbf{h}_{t+1}|\mathbf{h}_{t})-Kn\log(\Delta
t)$ (23)
$\displaystyle=\sum_{k=1}^{K}\Big{(}\log\mathrm{Pr}(h_{k,1})+\sum_{t=1}^{n-1}\log\mathrm{Pr}(h_{k,t+1}|h_{k,t})\Big{)}-Kn\log(\Delta
t)$ (24)

where (23)–(24) follow from the independencies described in (5).

Figure 10: Graphical model for (a) Polynomial and (b) GAN with added white
noise.

###  A.3 GAN

We approximate the GAN’s likelihood by calculating the likelihood of a model
which functions and gives results almost identical to the GAN (one with a
small amount of white noise added) using importance sampling and the
reparameterization method as in the Variational Autoencoder Kingma and Welling
(2013). Here, (9) is altered to

$\displaystyle X_{k,t+1}$
$\displaystyle=X_{k,t}+\omega_{k}(\mathbf{X}_{t})-\Delta t\
u(X_{k,t},h_{k,t+1})$ (25)

where $u(X_{k,t},h_{k,t+1})=U(X_{k,t},h_{k,t+1})+N_{k,t+1}$ and
$N_{k,t}\sim\mathcal{N}(0,0.001^{2})$. The log-likelihood is

$\displaystyle\log\mathrm{Pr}(\mathbf{x}_{1},...,\mathbf{x}_{n};\mathbf{x}_{0})$
$\displaystyle=\log\mathbb{E}_{\mathbf{h}\sim\mathrm{Pr}_{\mathbf{\tilde{H}}}.}\frac{\mathrm{Pr}(\mathbf{u}.|\mathbf{h}.,\mathbf{x}_{0})\mathrm{Pr}_{\mathbf{H.}}(\mathbf{h}.)}{\mathrm{Pr}_{\mathbf{\tilde{H}.}}(\mathbf{h}.)}-Kn\log(\Delta
t)$

where the expectation can now be approximated by a large sum if there is a
suitable importance sampler $\mathrm{Pr}_{\mathbf{\tilde{H}.}}(\mathbf{h}.)$.
50 samples were used for the importance sampling, and this was repeated 50
times to give 95% confidence intervals. We note that given the finite number
of samples taken it is of course possible that the GAN likelihood is larger,
but a well-trained importance sampler should minimise the chance of this. The
derivation and form of the importance sampler is detailed below.

The graphical model of the GAN with a small amount of white noise added is
shown in Figure 10b. $\mathbf{X}_{t}$ is related to $\mathbf{u}_{t}$ as in
(25) so

$\displaystyle(\mathbf{X}_{t}|\mathbf{X}_{t-1}{=}\mathbf{x}_{t-1},...,\mathbf{X}_{1}{=}\mathbf{x}_{1};\mathbf{X}_{0}{=}\mathbf{x}_{0})$
$\displaystyle=f(\mathbf{x}_{t-1})-\Delta
t(\mathbf{u}_{t}|\mathbf{X}_{t-1}{=}\mathbf{x}_{t-1},...,\mathbf{X}_{1}{=}\mathbf{x}_{1};\mathbf{X}_{0}{=}\mathbf{x}_{0})$
(26)

and a change of variables for likelihood gives

$\log\mathrm{Pr}(\mathbf{x}_{t}|\mathbf{x}_{t-1},...,\mathbf{x}_{1};\mathbf{x}_{0})=\log\mathrm{Pr}(\mathbf{u}_{t}|\mathbf{u}_{t-1},...,\mathbf{u}_{1};\mathbf{x}_{0})-K\log(\Delta
t)$ (27)

(13) is therefore

$\displaystyle\begin{split}&=\log\mathrm{Pr}(\mathbf{u}_{1}|\mathbf{x}_{0})+\log\mathrm{Pr}(\mathbf{u}_{2}|\mathbf{u}_{1},\mathbf{x}_{0})+\log\mathrm{Pr}(\mathbf{u}_{3}|\mathbf{u}_{2},\mathbf{u}_{1},\mathbf{x}_{0})\\\
{}&\qquad\qquad+...+\log\mathrm{Pr}(\mathbf{u}_{n}|\mathbf{u}_{n-1:0},\mathbf{x}_{0})-Kn\log(\Delta
t)\end{split}$
$\displaystyle=\log\mathrm{Pr}(\mathbf{u}.|\mathbf{x}_{0})-Kn\log(\Delta t)$

and just decomposing the first term below:

$\displaystyle\log\mathrm{Pr}(\mathbf{u}.|\mathbf{x}_{0})$
$\displaystyle=\log\mathbb{E}_{\mathbf{h}\sim\mathrm{Pr}_{\mathbf{H}}.}\mathrm{Pr}(\mathbf{u}.|\mathbf{h}.,\mathbf{x}_{0})$
$\displaystyle=\log\mathbb{E}_{\mathbf{h}\sim\mathrm{Pr}_{\mathbf{\tilde{H}}}.}\frac{\mathrm{Pr}(\mathbf{u}.|\mathbf{h}.,\mathbf{x}_{0})\mathrm{Pr}_{\mathbf{H.}}(\mathbf{h}.)}{\mathrm{Pr}_{\mathbf{\tilde{H}.}}(\mathbf{h}.)}$
(28)

where (28) can be approximated with a sufficiently large sum given a good
enough encoder $\mathrm{Pr}_{\mathbf{\tilde{H}.}}(\mathbf{h}.)$.

For training purposes, the lower-bound is used:

$\displaystyle\geq\mathbb{E}_{\mathbf{h}\sim\mathrm{Pr}_{\mathbf{\tilde{H}}}.}\log\frac{\mathrm{Pr}(\mathbf{u}.|\mathbf{h}.,\mathbf{x}_{0})\mathrm{Pr}_{\mathbf{H.}}(\mathbf{h}.)}{\mathrm{Pr}_{\mathbf{\tilde{H}.}}(\mathbf{h}.)}$
(29)

The terms in the numerator decompose as follows, using the independencies from
the graphical model and associated equations:

$\displaystyle\log\mathrm{Pr}(\mathbf{u}.|\mathbf{h}.,\mathbf{x}_{0})$
$\displaystyle=\sum_{t=1}^{n}\log\mathrm{Pr}(\mathbf{u}_{t}|\mathbf{h}_{t},\mathbf{x}_{t-1})$
$\displaystyle=\sum_{t=1}^{n}\sum_{k=1}^{K}\log\mathrm{Pr}(u_{k,t}|h_{k,t},x_{k,t-1})$

and

$\displaystyle\log\mathrm{Pr}_{\mathbf{H.}}(\mathbf{h}.)$
$\displaystyle=\log\mathrm{Pr}_{\mathbf{H}}(\mathbf{h}_{1})+\sum_{t=2}^{n}\log\mathrm{Pr}_{\mathbf{H}}(\mathbf{h}_{t}|\mathbf{h}_{t-1})$
$\displaystyle=\sum_{k=1}^{K}\Big{(}\log\mathrm{Pr}_{H}(h_{k,1})+\sum_{t=2}^{n}\log\mathrm{Pr}_{H}(h_{k,t}|h_{k,t-1})\Big{)}$

The perfect importance sampling distribution would be proportional to
$\mathrm{Pr}(\mathbf{h}.|\mathbf{u}.,\mathbf{x}_{0})$. Therefore we want
$\mathrm{Pr}_{\mathbf{\tilde{H}.}}(\mathbf{h}.)$ to allow for the same
dependencies. We choose it by carrying out the following steps:

$\displaystyle\log\mathrm{Pr}_{\mathbf{\tilde{H}.}}(\mathbf{h}.)$
$\displaystyle\propto\log\mathrm{Pr}(\mathbf{h}.|\mathbf{u}.,\mathbf{x}_{0})$
$\displaystyle=\log\mathrm{Pr}_{\mathbf{H}_{1}}(\mathbf{h}_{1}|\mathbf{u}.,\mathbf{x}_{0})+\sum_{t=2}^{n}\log\mathrm{Pr}_{\mathbf{H}_{t}}(\mathbf{h}_{t}|\mathbf{h}_{t-1},\mathbf{u}.,\mathbf{x}_{0})$
(30)

where (30) follows from using the laws of probability to decompose the term
above, followed by applications of the independencies from the graphical
model. Note that although the process is Markov it is not necessarily time-
homogeneous. To deal with this, an RNN is used to model
$\mathrm{Pr}_{\mathbf{H}_{t}}(\mathbf{h}_{t}|\mathbf{h}_{t-1},\mathbf{u}.,\mathbf{x}_{0})$
as
$\mathrm{Pr}_{\mathbf{H}}(\mathbf{h}_{t}|\mathbf{h}_{t-1},\mathbf{u}.,\mathbf{x}_{0},\mathbf{s}_{t})$
where the state $\mathbf{s}_{t}$ can in principle account for the non-
homogeneous updating procedure of $\mathbf{h}_{t}$. The log-likelihood is
therefore

$\displaystyle=\sum_{k=1}^{K}\Big{(}\log\mathrm{Pr}_{H_{1}}(h_{k,1}|u_{k.},x_{k,0})+\sum_{t=2}^{n}\log\mathrm{Pr}_{H_{t}}(h_{k,t}|h_{k,t-1},u_{k.},x_{k.},s_{k,t})\Big{)}$
(31)

where the spatial independencies introduced in (31) follow from the graphical
model and $u_{k,.}=u_{k,1:n}$ and $x_{k,.}=x_{k,0:n-1}$. The only simplifying
assumption used here is that the same update model for $h_{k,t}$ is used for
each component $k$. Finally, the sequence of $(u_{k.},x_{k.})$ is summarised
using a bidirectional RNN to give $w_{k,t}$. Therefore, the final form of the
log-likelihood of our importance sampler is

$\displaystyle=\sum_{k=1}^{K}\Big{(}\log\mathrm{Pr}_{H_{1}}(h_{k,1}|w_{k,1},x_{k,0})+\sum_{t=2}^{n}\log\mathrm{Pr}_{H_{t}}(h_{k,t}|h_{k,t-1},w_{k,t},x_{k,t-1},s_{k,t})\Big{)}$
(32)

The training of the importance sampler is done by maximising (29), with the
learnable parameters being those of the model used to learn $h_{k,1}$, the RNN
used to model the update of $h_{k,t}$, and the bidirectional RNN.

All authors supported the design of the study. RP created the models and
conducted research. RP prepared the manuscript with contributions from all co-
authors.

The authors declare that they have no conflict of interest.

###### Acknowledgements.

R.P. was funded by the Engineering and Physical Sciences Research Council
[grant number EP/S022961/1]. H.M.C. was funded by Natural Environment Research
Council [grant number NE/P018238/1]. J.S.H. was supported by the British
Antarctic Survey, Natural Environment Research Council (NERC) National
Capability funding. D.J.W. was funded by the University of Cambridge.

## References

  * Arcomano et al. (2022) Arcomano, T., Szunyogh, I., Wikner, A., Pathak, J., Hunt, B. R., and Ott, E.: A Hybrid Approach to Atmospheric Modeling That Combines Machine Learning With a Physics-Based Numerical Model, Journal of Advances in Modeling Earth Systems, 14, e2021MS002 712, 2022.
  * Arjovsky et al. (2017) Arjovsky, M., Chintala, S., and Bottou, L.: Wasserstein generative adversarial networks, in: International conference on machine learning, pp. 214–223, PMLR, 2017.
  * Arnold et al. (2013) Arnold, H. M., Moroz, I. M., and Palmer, T. N.: Stochastic parametrizations and model uncertainty in the Lorenz’96 system, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 371, 20110 479, 2013.
  * Berner et al. (2012) Berner, J., Jung, T., and Palmer, T.: Systematic model error: The impact of increased horizontal resolution versus improved stochastic and deterministic parameterizations, Journal of Climate, 25, 4946–4962, 2012.
  * Beucler et al. (2020) Beucler, T., Pritchard, M., Gentine, P., and Rasp, S.: Towards physically-consistent, data-driven models of convection, in: IGARSS 2020-2020 IEEE International Geoscience and Remote Sensing Symposium, pp. 3987–3990, IEEE, 2020.
  * Beucler et al. (2021) Beucler, T., Pritchard, M., Rasp, S., Ott, J., Baldi, P., and Gentine, P.: Enforcing analytic constraints in neural networks emulating physical systems, Physical Review Letters, 126, 098 302, 2021.
  * Bolton and Zanna (2019) Bolton, T. and Zanna, L.: Applications of deep learning to ocean data inference and subgrid parameterization, Journal of Advances in Modeling Earth Systems, 11, 376–399, 2019.
  * Brenowitz and Bretherton (2018) Brenowitz, N. D. and Bretherton, C. S.: Prognostic validation of a neural network unified physics parameterization, Geophysical Research Letters, 45, 6289–6298, 2018.
  * Brenowitz and Bretherton (2019) Brenowitz, N. D. and Bretherton, C. S.: Spatially extended tests of a neural network parametrization trained by coarse-graining, Journal of Advances in Modeling Earth Systems, 11, 2728–2744, 2019.
  * Buizza et al. (1999) Buizza, R., Milleer, M., and Palmer, T. N.: Stochastic representation of model uncertainties in the ECMWF ensemble prediction system, Quarterly Journal of the Royal Meteorological Society, 125, 2887–2908, 1999.
  * Chattopadhyay et al. (2020a) Chattopadhyay, A., Hassanzadeh, P., and Subramanian, D.: Data-driven predictions of a multiscale Lorenz 96 chaotic system using machine-learning methods: reservoir computing, artificial neural network, and long short-term memory network, Nonlinear Processes in Geophysics, 27, 373–389, 2020a.
  * Chattopadhyay et al. (2020b) Chattopadhyay, A., Subel, A., and Hassanzadeh, P.: Data-driven super-parameterization using deep learning: Experimentation with multiscale Lorenz 96 systems and transfer learning, Journal of Advances in Modeling Earth Systems, 12, e2020MS002 084, 2020b.
  * Cho et al. (2014) Cho, K., Van Merriënboer, B., Gulcehre, C., Bahdanau, D., Bougares, F., Schwenk, H., and Bengio, Y.: Learning phrase representations using RNN encoder-decoder for statistical machine translation, arXiv preprint arXiv:1406.1078, 2014.
  * Christensen et al. (2015) Christensen, H. M., Moroz, I. M., and Palmer, T. N.: Simulating weather regimes: Impact of stochastic and perturbed parameter schemes in a simple atmospheric model, Climate Dynamics, 44, 2195–2214, 2015.
  * Christensen et al. (2017) Christensen, H. M., Berner, J., Coleman, D. R., and Palmer, T.: Stochastic parameterization and El Niño–Southern Oscillation, Journal of Climate, 30, 17–38, 2017.
  * Chung et al. (2016) Chung, J., Ahn, S., and Bengio, Y.: Hierarchical multiscale recurrent neural networks, arXiv preprint arXiv:1609.01704, 2016.
  * Crommelin and Vanden-Eijnden (2008) Crommelin, D. and Vanden-Eijnden, E.: Subgrid-scale parameterization with conditional Markov chains, Journal of the Atmospheric Sciences, 65, 2661–2675, 2008.
  * Dawson and Palmer (2015) Dawson, A. and Palmer, T.: Simulating weather regimes: Impact of model resolution and stochastic parameterization, Climate Dynamics, 44, 2177–2193, 2015\.
  * Eck and Schmidhuber (2002) Eck, D. and Schmidhuber, J.: A first look at music composition using lstm recurrent neural networks, Istituto Dalle Molle Di Studi Sull Intelligenza Artificiale, 103, 48, 2002.
  * Esteban et al. (2017) Esteban, C., Hyland, S. L., and Rätsch, G.: Real-valued (medical) time series generation with recurrent conditional gans, arXiv preprint arXiv:1706.02633, 2017.
  * Gagne et al. (2020) Gagne, D. J., Christensen, H. M., Subramanian, A. C., and Monahan, A. H.: Machine learning for stochastic parameterization: Generative adversarial networks in the Lorenz’96 model, Journal of Advances in Modeling Earth Systems, 12, e2019MS001 896, 2020.
  * Gentine et al. (2018) Gentine, P., Pritchard, M., Rasp, S., Reinaudi, G., and Yacalis, G.: Could machine learning break the convection parameterization deadlock?, Geophysical Research Letters, 45, 5742–5751, 2018.
  * Goodfellow et al. (2014) Goodfellow, I. J., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Courville, A., and Bengio, Y.: Generative adversarial networks, arXiv preprint arXiv:1406.2661, 2014.
  * Graves (2013) Graves, A.: Generating sequences with recurrent neural networks, arXiv preprint arXiv:1308.0850, 2013.
  * Grover et al. (2018) Grover, A., Dhar, M., and Ermon, S.: Flow-gan: Combining maximum likelihood and adversarial learning in generative models, in: Proceedings of the AAAI Conference on Artificial Intelligence, vol. 32, 2018.
  * Guillaumin and Zanna (2021) Guillaumin, A. P. and Zanna, L.: Stochastic-deep learning parameterization of ocean momentum forcing, Journal of Advances in Modeling Earth Systems, 13, e2021MS002 534, 2021.
  * Gulrajani et al. (2017) Gulrajani, I., Ahmed, F., Arjovsky, M., Dumoulin, V., and Courville, A.: Improved training of wasserstein gans, arXiv preprint arXiv:1704.00028, 2017.
  * Hasselmann (1976) Hasselmann, K.: Stochastic climate models part I. Theory, tellus, 28, 473–485, 1976\.
  * Hochreiter and Schmidhuber (1997) Hochreiter, S. and Schmidhuber, J.: Long short-term memory, Neural computation, 9, 1735–1780, 1997.
  * Johnson et al. (2019) Johnson, S. J., Stockdale, T. N., Ferranti, L., Balmaseda, M. A., Molteni, F., Magnusson, L., Tietsche, S., Decremer, D., Weisheimer, A., Balsamo, G., et al.: SEAS5: the new ECMWF seasonal forecast system, Geoscientific Model Development, 12, 1087–1117, 2019.
  * Kingma and Ba (2014) Kingma, D. P. and Ba, J.: Adam: A method for stochastic optimization, arXiv preprint arXiv:1412.6980, 2014.
  * Kingma and Welling (2013) Kingma, D. P. and Welling, M.: Auto-encoding variational bayes, arXiv preprint arXiv:1312.6114, 2013.
  * Krasnopolsky et al. (2013) Krasnopolsky, V. M., Fox-Rabinovitz, M. S., and Belochitski, A. A.: Using ensemble of neural networks to learn stochastic convection parameterizations for climate and numerical weather prediction models from data simulated by a cloud resolving model, Advances in Artificial Neural Systems, 2013, 2013.
  * Kwasniok (2012) Kwasniok, F.: Data-based stochastic subgrid-scale parametrization: an approach using cluster-weighted modelling, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 370, 1061–1086, 2012\.
  * Leutbecher and Palmer (2008) Leutbecher, M. and Palmer, T. N.: Ensemble forecasting, Journal of computational physics, 227, 3515–3539, 2008.
  * Leutbecher et al. (2017) Leutbecher, M., Lock, S.-J., Ollinaho, P., Lang, S. T., Balsamo, G., Bechtold, P., Bonavita, M., Christensen, H. M., Diamantakis, M., Dutra, E., et al.: Stochastic representations of model uncertainties at ECMWF: State of the art and future vision, Quarterly Journal of the Royal Meteorological Society, 143, 2315–2339, 2017.
  * Liu et al. (2020) Liu, Y., Kutz, J. N., and Brunton, S. L.: Hierarchical deep learning of multiscale differential equation time-steppers, arXiv preprint arXiv:2008.09768, 2020.
  * Lorenz (1996) Lorenz, E. N.: Predictability: A problem partly solved, in: Proc. Seminar on predictability, vol. 1, 1996.
  * Lorenz (2006) Lorenz, E. N.: Regimes in simple systems, Journal of the atmospheric sciences, 63, 2056–2073, 2006.
  * Mogren (2016) Mogren, O.: C-RNN-GAN: Continuous recurrent neural networks with adversarial training, arXiv preprint arXiv:1611.09904, 2016.
  * Molteni et al. (2011) Molteni, F., Stockdale, T., Balmaseda, M., Balsamo, G., Buizza, R., Ferranti, L., Magnusson, L., Mogensen, K., Palmer, T., and Vitart, F.: The new ECMWF seasonal forecast system (System 4), vol. 49, European Centre for medium-range weather forecasts Reading, 2011.
  * O’Gorman and Dwyer (2018) O’Gorman, P. A. and Dwyer, J. G.: Using machine learning to parameterize moist convection: Potential for modeling of climate, climate change, and extreme events, Journal of Advances in Modeling Earth Systems, 10, 2548–2563, 2018.
  * Palmer (2012) Palmer, T. N.: Towards the probabilistic Earth-system simulator: A vision for the future of climate and weather prediction, Quarterly Journal of the Royal Meteorological Society, 138, 841–861, 2012.
  * Palmer et al. (2009) Palmer, T. N., Buizza, R., Doblas-Reyes, F., Jung, T., Leutbecher, M., Shutts, G. J., Steinheimer, M., and Weisheimer, A.: Stochastic parametrization and model uncertainty, 2009.
  * Rasp (2020) Rasp, S.: Coupled online learning as a way to tackle instabilities and biases in neural network parameterizations: general algorithms and Lorenz 96 case study (v1. 0), Geoscientific Model Development, 13, 2185–2196, 2020.
  * Rasp et al. (2018) Rasp, S., Pritchard, M. S., and Gentine, P.: Deep learning to represent subgrid processes in climate models, Proceedings of the National Academy of Sciences, 115, 9684–9689, 2018.
  * Sanchez et al. (2016) Sanchez, C., Williams, K. D., and Collins, M.: Improved stochastic physics schemes for global weather and climate models, Quarterly Journal of the Royal Meteorological Society, 142, 147–159, 2016.
  * Schneider et al. (2017) Schneider, T., Teixeira, J., Bretherton, C. S., Brient, F., Pressel, K. G., Schär, C., and Siebesma, A. P.: Climate goals and computing the future of clouds, Nature Climate Change, 7, 3–5, 2017.
  * Selten and Branstator (2004) Selten, F. M. and Branstator, G.: Preferred regime transition routes and evidence for an unstable periodic orbit in a baroclinic model, Journal of the atmospheric sciences, 61, 2267–2282, 2004.
  * Skamarock et al. (2019) Skamarock, W. C., Klemp, J. B., Dudhia, J., Gill, D. O., Liu, Z., Berner, J., Wang, W., Powers, J. G., Duda, M. G., Barker, D. M., et al.: A description of the advanced research WRF model version 4, National Center for Atmospheric Research: Boulder, CO, USA, 145, 145, 2019.
  * Stephenson et al. (2004) Stephenson, D., Hannachi, A., and O’Neill, A.: On the existence of multiple climate regimes, Quarterly Journal of the Royal Meteorological Society: A journal of the atmospheric sciences, applied meteorology and physical oceanography, 130, 583–605, 2004.
  * Stockdale et al. (2011) Stockdale, T. N., Anderson, D. L., Balmaseda, M. A., Doblas-Reyes, F., Ferranti, L., Mogensen, K., Palmer, T. N., Molteni, F., and Vitart, F.: ECMWF seasonal forecast system 3 and its prediction of sea surface temperature, Climate dynamics, 37, 455–471, 2011.
  * Straus et al. (2007) Straus, D. M., Corti, S., and Molteni, F.: Circulation regimes: Chaotic variability versus SST-forced predictability, Journal of climate, 20, 2251–2272, 2007.
  * Sutskever et al. (2011) Sutskever, I., Martens, J., and Hinton, G. E.: Generating text with recurrent neural networks, in: ICML, 2011.
  * Sutskever et al. (2014) Sutskever, I., Vinyals, O., and Le, Q. V.: Sequence to sequence learning with neural networks, in: Advances in neural information processing systems, pp. 3104–3112, 2014.
  * Theis et al. (2015) Theis, L., Oord, A. v. d., and Bethge, M.: A note on the evaluation of generative models, arXiv preprint arXiv:1511.01844, 2015.
  * Vaswani et al. (2017) Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, Ł., and Polosukhin, I.: Attention is all you need, in: Advances in neural information processing systems, pp. 5998–6008, 2017.
  * Vlachas et al. (2018) Vlachas, P. R., Byeon, W., Wan, Z. Y., Sapsis, T. P., and Koumoutsakos, P.: Data-driven forecasting of high-dimensional chaotic systems with long short-term memory networks, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 474, 20170 844, 2018.
  * Vlachas et al. (2020) Vlachas, P. R., Pathak, J., Hunt, B. R., Sapsis, T. P., Girvan, M., Ott, E., and Koumoutsakos, P.: Backpropagation algorithms and reservoir computing in recurrent neural networks for the forecasting of complex spatiotemporal dynamics, Neural Networks, 126, 191–217, 2020.
  * Vlachas et al. (2022) Vlachas, P. R., Arampatzis, G., Uhler, C., and Koumoutsakos, P.: Multiscale simulations of complex systems by learning their effective dynamics, Nature Machine Intelligence, 4, 359–366, 2022.
  * Walters et al. (2019) Walters, D., Baran, A. J., Boutle, I., Brooks, M., Earnshaw, P., Edwards, J., Furtado, K., Hill, P., Lock, A., Manners, J., et al.: The Met Office Unified Model global atmosphere 7.0/7.1 and JULES global land 7.0 configurations, Geoscientific Model Development, 12, 1909–1963, 2019.
  * Yuval and O’Gorman (2020) Yuval, J. and O’Gorman, P. A.: Stable machine-learning parameterization of subgrid processes for climate modeling at a range of resolutions, Nature communications, 11, 1–10, 2020.
  * Yuval et al. (2021) Yuval, J., O’Gorman, P. A., and Hill, C. N.: Use of neural networks for stable, accurate and physically consistent parameterization of subgrid atmospheric processes with good performance at reduced precision, Geophysical Research Letters, 48, e2020GL091 363, 2021.
  * Zelinka et al. (2020) Zelinka, M. D., Myers, T. A., McCoy, D. T., Po-Chedley, S., Caldwell, P. M., Ceppi, P., Klein, S. A., and Taylor, K. E.: Causes of higher climate sensitivity in CMIP6 models, Geophysical Research Letters, 47, e2019GL085 782, 2020.
  * Zhao et al. (2020) Zhao, M., Cong, Y., Dai, S., and Carin, L.: Bridging Maximum Likelihood and Adversarial Learning via $\alpha$-Divergence, in: Proceedings of the AAAI Conference on Artificial Intelligence, vol. 34, pp. 6901–6908, 2020.