# Invalid Logic, Equivalent Gains:
The Bizarreness of Reasoning in Language Model Prompting

Rylan Schaeffer  Kateryna Pistunova  Samar Khanna  Sarthak Consul  Sanmi
Koyejo

###### Abstract

Language models can be prompted to reason through problems in a manner that
significantly improves performance. However, why such prompting improves
performance is unclear. Recent work showed that using logically invalid Chain-
of-Thought (CoT) prompting improves performance almost as much as logically
valid CoT prompting, and that editing CoT prompts to replace problem-specific
information with abstract information or out-of-distribution information
typically doesn’t harm performance. Critics have responded that these findings
are based on too few and too easily solved tasks to draw meaningful
conclusions. To resolve this dispute, we test whether logically invalid CoT
prompts offer the same level of performance gains as logically valid prompts
on the hardest tasks in the BIG-Bench benchmark, termed BIG-Bench Hard (BBH).
We find that the logically invalid reasoning prompts do indeed achieve similar
performance gains on BBH tasks as logically valid reasoning prompts. We also
discover that some CoT prompts used by previous works contain logical errors.
This suggests that covariates beyond logically valid reasoning are responsible
for performance improvements.

Reasoning, Large Language Models, Language Models, ICML, Prompting, Chain-of-
Thought

##  1 Introduction

Language models can perform significantly better when prompted in particular
ways. For example, prompts that recommend or guide language models through
step-by-step processing have been shown to significantly improve performance
on question answering, conversational response generation and other tasks (Nye
et al., 2021; Wei et al., 2022b; Jung et al., 2022; Kojima et al., 2022; Yao
et al., 2023). These prompting techniques are especially effective on the
hardest tasks (Suzgun et al., 2022) in the BIG-Bench benchmark (Srivastava et
al., 2022), leading many to conclude that such techniques unlock
emergent111Although see Schaeffer et al. (2023). human-like reasoning
abilities in large language models (Wei et al., 2022a). However, why such
prompting strategies improve performance is unclear. Madaan & Yazdanbakhsh
(2022) showed replacing problem-specific information in Chain-of-Thought (CoT)
prompts with either abstract information or out-of-distribution information
typically doesn’t harm CoT’s performance gains, and Wang et al. (2022) showed
that logically invalid CoT prompts (i.e., prompts which contain logically
invalid chains of reasoning) achieve approximately the same gains as logically
valid CoT prompts. These discoveries are puzzling, but critics warned against
drawing confident or far-ranging conclusions results because (1) only a small
number of tasks are considered (specifically GSM8K (Cobbe et al., 2021),
Bamboogle (Press et al., 2022), and two commonsense tasks (Srivastava et al.,
2022)) and (2) the tasks are relatively easy to solve.

In this work, we aim to clarify this dispute by asking whether logically
invalid CoT prompts indeed offer comparable performance gains as logically
valid CoT prompts in more and harder tasks. To answer this question, we
evaluate the performance of logically invalid CoT prompts on the hardest tasks
in BIG-Bench (Srivastava et al., 2022), the so-called BIG-Bench Hard (BBH)
tasks (Suzgun et al., 2022). We find that logically invalid CoT prompts induce
performance gains approaching logically valid CoT prompts on BBH tasks. We
also discover that CoT prompts from previous work contain logical errors, a
discovery that we confirmed with the original authors. Our findings support
growing evidence that while prompting techniques can, without doubt,
significantly improve language model performance on complex tasks, there is
questionable evidence that these performance gains are attributable to logical
reasoning skills. Covariates in the prompts other than logical reasoning may
be responsible for the performance improvements and point to the need for
continued investigation.

##  2 Methods

###  2.1 Tasks: BIG-Bench Hard $\subset$ BIG-Bench

Beyond-the-Imitation Game Benchmark (BIG-Bench) is a benchmark of over 200
natural language tasks created to evaluate the capabilities and limitations of
language models (Srivastava et al., 2022). BIG-Bench Hard (BBH) is a subset of
BIG-Bench consisting of 23 tasks specifically identified as extremely
difficult for language models (Suzgun et al., 2022). BBH was constructed by
excluding BIG-Bench tasks using the following criteria: tasks with more than
three subtasks, tasks with fewer than 103 examples, tasks without human-rater
baselines, tasks that do not use multiple-choice or exact match as the
evaluation metric, and tasks on which the best reported model beats average
reported human-rater score; the remaining 23 BIG-Bench tasks comprise the BBH
tasks. The BBH tasks cover a variety of categories, including traditional NLP,
mathematics, commonsense reasoning, and question-answering. In general, BBH
tasks typically fall into one of two categories: more traditional NLP tasks
(e.g., Date Understanding) or more algorithmic tasks (e.g., Boolean
Expressions). Many state of the art language models, including GPT-3 (Brown et
al., 2020) and PaLM (Chowdhery et al., 2022), as well as internal dense and
sparse Google models, were shown to be incapable of solving BBH tasks above
average human rater performance when asked directly.

###  2.2 Prompt Types

In this work, we compare three prompt types applied to the BBH tasks. Examples
of task problems and prompts are shown in Table 1.

#### Chain-of-Thought (CoT) Prompting

In (logically valid) Chain-of-Thought (CoT) prompting, each question in each
task is prepended with three high-quality, human hand-written, topical
examples of Question$\rightarrow$Reasoning$\rightarrow$Answer. The reasoning
is a sequence of natural language steps that derive the answer from
information in the question (Nye et al., 2021; Wei et al., 2022b). For each
BBH task, Suzgun et al. (2022) released three examples of
Question$\rightarrow$Reasoning$\rightarrow$Answer that we adapt. We do not use
their prompts directly because we want CoT prompts to contain logically valid
chains of reasoning, and as we discovered (Sec. 3.1), some of Suzgun et al.’s
CoT prompts contain logical errors.

#### Logically Invalid Chain-of-Thought (Invalid CoT) Prompting

In logically invalid Chain-of-Thought (Invalid CoT) prompting, we edit the
reasoning in the three examples for all three tasks to become logically
invalid. To do this, we modify the reasoning to contain nonsensical or invalid
steps. For instance, on the Boolean Expressions task, we provided the
reasoning, “Because English does not permit multiple negatives, the expression
”(not not True)” evaluates to ”( not True )”, and on the Date Understanding
task, we provided the reasoning, “If tomorrow is 11/12/2019, then today is
11/12/2018. The date one year ago from today is 11/11/2018.” In all cases, the
modified BBH exemplars were logically invalid, but reached the correct
answers. We wrote each example once and did not modify the examples whatsoever
to improve task performance.

#### Answer-Only (AO) Prompting

To establish baseline performance, we follow previous work and prompt language
models to answer each question in each task directly. We call this “Answer
Only” (AO) prompting.

Prompt Type |  Example Query

Evaluate the result of a random Boolean expression.

Q: not ( ( not not True ) ) is
---|---
Answer Only |  False
Scratchpad/Chain-of-Thought (CoT) |  A: Let’s think step by step. Remember that (i) expressions inside brackets are always evaluated first and that (ii) the order of operations from highest priority to lowest priority is ”not”, ”and”, ”or”, respectively. We first simplify this expression ”Z” as follows: ”Z = not ( ( not not True ) ) = not ( ( A ) )” where ”A = not not True”. Let’s evaluate A: A = not not True = not (not True) = not False = True. Plugging in A, we get: Z = not ( ( A ) ) = not ( ( True ) ) = not True = False. So the answer is False.
Logically-Invalid Chain-of-Thought (Invalid CoT) |  A: Let’s think step by step. Remember that (i) expressions inside brackets are always evaluated first and that (ii) the order of operations from highest priority to lowest priority is ”not”, ”and”, ”or”, respectively. Because English does not permit multiple negatives, the expression ”(not not True)” evaluates to ”( not True )”. The expression ”not ( ( not not True ) )” therefore evaluates to ”not ( ( not True ) )”. By the same logic, the expression ”not ( ( not True ) )” simplifies to ”not True”. In Boolean logic, ”not True” is False. So the answer is False.
Table 1: Examples of different prompt types.

###  2.3 Metrics

To evaluate the performance of different prompting types, we used Accuracy
(i.e. Exact String Match) to evaluate model outputs. We acknowledge that
recent work showed such a choice of metric may reach misleading scientific
conclusions (Schaeffer et al., 2023), but made the choice to ensure fair
comparison with the relevant preceding work.

###  2.4 Choice of Language Model

Figure 1: Codex benefits more from Chain-of-Thought Prompting than InstructGPT
or PaLM 540B on BIG-Bench Hard Tasks. Dashed vertical lines indicate average
model performance across all BIG-Bench Hard Tasks. Data from Suzgun et al.
(2022).

Due to limited resources, we could only evaluate one model. Only Codex (Chen
et al., 2021) and InstructGPT (Ouyang et al., 2022) were publicly queryable,
thus ruling out PaLM 540B (Chowdhery et al., 2022). We chose to evaluate Codex
because our visualizing Suzgun et al’s data revealed Codex outperformed both
InstructGPT and PaLM 540B; this point was not made by the authors, but
visualizing their data anew reveals this conclusion (Fig. 1). Although this
finding may seem surprising, independent work has found that language models
reason better if pretrained heavily on code (Madaan et al., 2022), as Codex
was. Perhaps unsurprisingly, we found that Codex’s high average performance
with CoT prompting is driven by large improvements on BBH algorithmic tasks,
overriding mediocre improvements (or regressions) on natural language tasks.

###  2.5 Choice of BBH Tasks

We evaluate AO, CoT, and Invalid CoT prompting strategies on all BBH tasks.

##  3 Experimental Results

Figure 2: Sanity checks. Top: All prompting strategies almost always produced
an answer across tasks. Bottom: Our accuracies for Answer Only (AO) and Chain-
of-Thought (CoT) prompting closely matched published values from Suzgun et al.
(2022) except on one task: word_sorting.

###  3.1 Previous CoT Prompts Contain Logical Errors

While converting the logically valid CoT prompts written for BBH tasks by
Suzgun et al. (2022) to logically invalid CoT prompts, we discovered that at
least three of the BBH tasks’ CoT prompts were already logically invalid. On
the Multistep Arithmetic Task, a human handwritten CoT prompt reads: “Then,
the final equation is A * B = -41 * -3 = (-61) * (-3) = 123. So the answer is
123.” To clarify, -41 was substituted with -61 without justification, and -61
* -3 is not 123, even though 123 is the correct answer. On the Navigate task,
a human handwritten CoT prompt reads: “(4) Take 7 steps right: (0, 7), facing
the positive y-axis.” should be “(4) Take 7 steps right: (0, 0), facing the
positive y-axis.” On the Web of Lies task, “(3) Vina says Jerry tells the
truth. Since we know from (2) that Jerry tells the truth, if Vine says Jerry
tells the truth, then Vina tells the truth.” should be “(3) Vina says Jerry
tells the truth. Since we know from (2) that Jerry tells the truth, if Vina
says Jerry tells the truth, then Vina tells the truth.” All three errors were
raised to the authors, who confirmed our discoveries and accepted our
suggested corrections as well. This discovery already hints that CoT prompts
need not be logically correct for the model to output the correct answer. It
also hints that previous CoT prompts from earlier papers, e.g., (Wei et al.,
2022b), may need to be investigated further.

###  3.2 Sanity Checks

We found that Codex produced an answer under all prompting strategies on
almost all tasks (Fig. 2 top), and the accuracies under Answer Only (AO) and
Chain-of-Thought (CoT) prompting closely matched published values (Suzgun et
al., 2022). The one exception was a single task: word_sorting.

Figure 3: Accuracy per task per prompt type on BIG-Bench Hard. Prompt Types:
AO = Answer Only, CoT = Chain-of-Thought, CoT (Invalid) = Logically Invalid
Chain-of-Thought. Approximately 200-250 questions per BIG-Bench Hard task.

###  3.3 Accuracy by Prompt Type by Task

The Accuracy of Codex on each of the BIG-Bench Hard (BBH) tasks under each of
the three prompt types (Answer Only, Chain-of-Thought, Logically Invalid
Chain-of-Thought) is displayed in Fig. 3. Despite using the same model and the
same decoding hyperparameters (e.g., temperature), we notice that there is
tremendous variation in the accuracies. To better compare the performance
across the three prompt types, we visualized each prompt type’s average
accuracy over all BBH tasks. We found that Chain-of-Thought prompting beats
Answer Only prompting, but Logically Invalid Chain-of-Thought is close behind
Chain-of-Thought and better than Answer Only (Fig. 4 top). Pairwise
comparisons between prompt types across different BBH tasks reveal that
different prompt types do not dominate one another (Fig. 4 bottom), but rather
display rich structure. We note that we do not optimize logically invalid CoT
prompts whatsoever, meaning the Invalid CoT performance here is likely below
what tuning could achieve.

Figure 4: Logically invalid Chain-of-Thought outperforms Answer Only and
nearly matches Chain-of-Thought. Top: Accuracy per BBH task, averaged over all
BBH tasks. Bottom: Pairwise accuracy plots comparing the three different
prompt types.

##  4 Conclusion

On the diverse and challenging BIG-Bench Hard tasks, we find that Chain-of-
Thought prompting performs best on average, but logically invalid Chain-of-
Thought prompting is close behind and outperforms Answer Only prompting. This
demonstrates that completely illogical reasoning in the CoT prompts do not
significantly harm the performance of the language model. Our findings suggest
that valid reasoning in prompting is not the chief driver of performance
gains, raising the question of what is. We note that there are complementary
approaches towards achieving reasoning in language models such as enforcing
valid reasoning in the model architecture (Creswell et al., 2022; Creswell &
Shanahan, 2022) or via autoformalization (Wu et al., 2022; Azerbayev et al.,
2023).

Our work raises important questions for future work. Why are models robust to
invalid CoT prompts? What features of the data or prompts result in the model
outputting inconsistent or invalid outputs? Does increasing the degree of
”incorrectness” or the number of incorrect prompts affect the model’s
sensitivity to invalid CoT? What other properties of the valid prompts is the
model sensitive to? Answering these questions can yield useful insights into
prompt engineering for LLMs, as well as a deeper understanding of when models
output inconsistent or “hallucinated” answers.

## References

  * Azerbayev et al. (2023) Azerbayev, Z., Piotrowski, B., Schoelkopf, H., Ayers, E. W., Radev, D., and Avigad, J.  Proofnet: Autoformalizing and formally proving undergraduate-level mathematics.  _arXiv preprint arXiv:2302.12433_ , 2023.
  * Brown et al. (2020) Brown, T., Mann, B., Ryder, N., Subbiah, M., Kaplan, J. D., Dhariwal, P., Neelakantan, A., Shyam, P., Sastry, G., Askell, A., et al.  Language models are few-shot learners.  _Advances in neural information processing systems_ , 33:1877–1901, 2020.
  * Chen et al. (2021) Chen, M., Tworek, J., Jun, H., Yuan, Q., Pinto, H. P. d. O., Kaplan, J., Edwards, H., Burda, Y., Joseph, N., Brockman, G., et al.  Evaluating large language models trained on code.  _arXiv preprint arXiv:2107.03374_ , 2021.
  * Chowdhery et al. (2022) Chowdhery, A., Narang, S., Devlin, J., Bosma, M., Mishra, G., Roberts, A., Barham, P., Chung, H. W., Sutton, C., Gehrmann, S., et al.  Palm: Scaling language modeling with pathways.  _arXiv preprint arXiv:2204.02311_ , 2022.
  * Cobbe et al. (2021) Cobbe, K., Kosaraju, V., Bavarian, M., Chen, M., Jun, H., Kaiser, L., Plappert, M., Tworek, J., Hilton, J., Nakano, R., et al.  Training verifiers to solve math word problems.  _arXiv preprint arXiv:2110.14168_ , 2021.
  * Creswell & Shanahan (2022) Creswell, A. and Shanahan, M.  Faithful reasoning using large language models.  _arXiv preprint arXiv:2208.14271_ , 2022.
  * Creswell et al. (2022) Creswell, A., Shanahan, M., and Higgins, I.  Selection-inference: Exploiting large language models for interpretable logical reasoning.  _arXiv preprint arXiv:2205.09712_ , 2022.
  * Jung et al. (2022) Jung, J., Qin, L., Welleck, S., Brahman, F., Bhagavatula, C., Bras, R. L., and Choi, Y.  Maieutic prompting: Logically consistent reasoning with recursive explanations.  _arXiv preprint arXiv:2205.11822_ , 2022.
  * Kojima et al. (2022) Kojima, T., Gu, S. S., Reid, M., Matsuo, Y., and Iwasawa, Y.  Large language models are zero-shot reasoners.  _arXiv preprint arXiv:2205.11916_ , 2022.
  * Madaan & Yazdanbakhsh (2022) Madaan, A. and Yazdanbakhsh, A.  Text and patterns: For effective chain of thought, it takes two to tango.  _arXiv preprint arXiv:2209.07686_ , 2022.
  * Madaan et al. (2022) Madaan, A., Zhou, S., Alon, U., Yang, Y., and Neubig, G.  Language models of code are few-shot commonsense learners.  _arXiv preprint arXiv:2210.07128_ , 2022.
  * Nye et al. (2021) Nye, M., Andreassen, A. J., Gur-Ari, G., Michalewski, H., Austin, J., Bieber, D., Dohan, D., Lewkowycz, A., Bosma, M., Luan, D., et al.  Show your work: Scratchpads for intermediate computation with language models.  _arXiv preprint arXiv:2112.00114_ , 2021.
  * Ouyang et al. (2022) Ouyang, L., Wu, J., Jiang, X., Almeida, D., Wainwright, C., Mishkin, P., Zhang, C., Agarwal, S., Slama, K., Ray, A., et al.  Training language models to follow instructions with human feedback.  _Advances in Neural Information Processing Systems_ , 35:27730–27744, 2022.
  * Press et al. (2022) Press, O., Zhang, M., Min, S., Schmidt, L., Smith, N. A., and Lewis, M.  Measuring and narrowing the compositionality gap in language models.  _arXiv preprint arXiv:2210.03350_ , 2022.
  * Schaeffer et al. (2023) Schaeffer, R., Miranda, B., and Koyejo, S.  Are emergent abilities of large language models a mirage?  _arXiv preprint arXiv:2304.15004_ , 2023.
  * Srivastava et al. (2022) Srivastava, A., Rastogi, A., Rao, A., Shoeb, A. A. M., Abid, A., Fisch, A., Brown, A. R., Santoro, A., Gupta, A., Garriga-Alonso, A., et al.  Beyond the imitation game: Quantifying and extrapolating the capabilities of language models.  _arXiv preprint arXiv:2206.04615_ , 2022.
  * Suzgun et al. (2022) Suzgun, M., Scales, N., Schärli, N., Gehrmann, S., Tay, Y., Chung, H. W., Chowdhery, A., Le, Q. V., Chi, E. H., Zhou, D., et al.  Challenging big-bench tasks and whether chain-of-thought can solve them.  _arXiv preprint arXiv:2210.09261_ , 2022.
  * Wang et al. (2022) Wang, B., Min, S., Deng, X., Shen, J., Wu, Y., Zettlemoyer, L., and Sun, H.  Towards understanding chain-of-thought prompting: An empirical study of what matters.  _arXiv preprint arXiv:2212.10001_ , 2022.
  * Wei et al. (2022a) Wei, J., Tay, Y., Bommasani, R., Raffel, C., Zoph, B., Borgeaud, S., Yogatama, D., Bosma, M., Zhou, D., Metzler, D., et al.  Emergent abilities of large language models.  _arXiv preprint arXiv:2206.07682_ , 2022a.
  * Wei et al. (2022b) Wei, J., Wang, X., Schuurmans, D., Bosma, M., Chi, E., Le, Q., and Zhou, D.  Chain of thought prompting elicits reasoning in large language models.  _arXiv preprint arXiv:2201.11903_ , 2022b.
  * Wu et al. (2022) Wu, Y., Jiang, A. Q., Li, W., Rabe, M. N., Staats, C., Jamnik, M., and Szegedy, C.  Autoformalization with large language models.  _arXiv preprint arXiv:2205.12615_ , 2022.
  * Yao et al. (2023) Yao, S., Yu, D., Zhao, J., Shafran, I., Griffiths, T. L., Cao, Y., and Narasimhan, K.  Tree of thoughts: Deliberate problem solving with large language models.  _arXiv preprint arXiv:2305.10601_ , 2023.


[1]Xiaoxin Yin

[1]Independent Researcher, San Jose, California, USA

# Automated Discovery of Integral with Deep Learning

<EMAIL_ADDRESS>*

###### Abstract

Recent advancements in the realm of deep learning, particularly in the
development of large language models (LLMs), have demonstrated AI’s ability to
tackle complex mathematical problems or solving programming challenges.
However, the capability to solve well-defined problems based on extensive
training data differs significantly from the nuanced process of making
scientific discoveries. Trained on almost all human knowledge available,
today’s sophisticated LLMs basically learn to predict sequences of tokens.
They generate mathematical derivations and write code in a similar way as
writing an essay, and do not have the ability to pioneer scientific
discoveries in the manner a human scientist would do.

In this study we delve into the potential of using deep learning to rediscover
a fundamental mathematical concept: integrals. By defining integrals as area
under the curve, we illustrate how AI can deduce the integral of a given
function, exemplified by inferring $\int_{0}^{x}t^{2}dt=\frac{x^{3}}{3}$ and
$\int_{0}^{x}ae^{bt}dt=\frac{a}{b}e^{bx}-\frac{a}{b}$. Our experiments show
that deep learning models can approach the task of inferring integrals either
through a sequence-to-sequence model, akin to language translation, or by
uncovering the rudimentary principles of integration, such as
$\int_{0}^{x}t^{n}dt=\frac{x^{n+1}}{n+1}$.

###### keywords:

Deep Learning, Integral, Symbolic Regression

##  1 Introduction

The recent advances in deep learning, especially those in large language
models (aka LLM), have shown the possibility of an AI agent performing any
task a human can perform, including scientific research. Recent studies have
shown that LLMs such as GPT-4[1], Microsoft Copilot[2], and CodeLlama[3] can
solve competition-level coding problems [4], and LLMs such as GPT-4 and
Llemma[5] can solve some high-school-level competition math problems
(including some IMO-level problems). These LLMs can certainly help researchers
solve some problems they encounter in their daily research.

However, being able to solve a type of well-defined problems is very different
from making discoveries in scientific research. For instance, in order to
train an LLM to solve coding problems, a general-purpose LLM is often fine-
tuned on all public code on GitHub, and also fine-tuned on hundreds of
thousands of coding problems from various platforms such as CodeForce and
LeetCode. For example, CodeLlama-Python underwent fine-tuning with 100 billion
tokens of Python code. The LLM simply learns how to write code given the
coding problem (which is the prompt), by learning to predict the next token in
its code given the prompt and tokens it has generated. This is essentially the
same methodology used to train a model to write novels after reading millions
of novels. It does not have the capability of discovering what is not already
known to it, making it unable to make scientific discoveries like a scientist
would do.

Very recently DeepMind published their pioneer work on training a deep
learning model to solve IMO level geometry problems, without using theorems
found by human [6]. The authors randomly generated geometric data and let the
model learn the geometric rules from one billion examples. This is an exciting
piece of work in which an AI discovers mathematical theorems and methodologies
on its own.

In this paper we would like to explore how to use deep learning models to
discover an old but most important tool: Integrals. We define integral of a
function $f$ using area under the curve from 0 to $x$ (i.e., the integral of a
function $f(x)$ is defined as $\int_{0}^{x}f(t)dt$), which can be easily
calculated by sampling a function (e.g., using the simple Trapezoidal rule).
We hope to answer the following question: Simply given this definition, can
deep learning models and algorithms learn to derive the mathematical
expression of the integral of a function, and find the basic rules of
integrals for each type of functions (polynomials, exponential, trigonometry,
etc.)

In this study the AI models are only given the following data and tools:

  * •

A large set of randomly generated univariate functions.

  * •

The ability to evaluate a function.

  * •

The definition of integral (defined as area under the curve), and a program
that can calculate the area under the curve of any function using the
Trapezoidal rule, which mainly calculates the cumulative sum of a function.

Our experiment shows that given the above data and tools, AI models can
successfully finish the following two tasks:

  * •

Based on the mathematical expression of a function $f(x)$, the model can
directly infer the mathematical expression of the integral of $f(x)$, without
evaluating the function. The process is actually similar to translating a
sentence from one language to another one.

  * •

Discover the basic rules for integrating each simple type of functions, such
as polynomials, trigonometric, and exponential functions.

This experiment shows that AI models can acquire the capability of
integration, without human providing any examples. The only information
required is how integral is defined (as area under the curve in our study),
and AI can derive the rest. Please note that we train our AI models from
scratch and it never sees any information other than what we provide in our
experiment.

Our experiment is consisted of four steps. First we compute the definite
integral ($\int_{0}^{x}f(t)dt$) of each random function $f$. Then we use a
symbolic regression model[7] to infer the mathematical expression of the
definite integral $g(x)=\int_{0}^{x}f(t)dt$. Please note the symbolic
regression model was trained only on randomly generated functions and their
curves, and was not aware of the concept of integral. The third step is to
train a language model to infer the integral function $g$ from $f$ without
calculating area under the curve, just like what a college student would do in
the calculus class. Lastly, we try to infer rules of integral for each simple
type of functions. For example,
$\int_{0}^{x}ae^{bt}dt=\frac{a}{b}e^{bx}-\frac{a}{b}$, where $a$ and $b$ are
constants.

In this process we do not utilize any LLMs trained on human knowledge, because
they probably have learned about integrals during training. The only training
data we use are randomly generated functions, their curves, and their areas
under the curves.

We have no doubt that an AI agent would be more powerful if it has been
trained on human knowledge and thus is capable to perform reasoning like a
human. But that also means the agent has read all human knowledge, and we can
only test its discovery power on what humans have not discovered yet, which is
certainly a highly challenging problem. Therefore, in this study we choose to
study how AI can discover an important math concept like integral, without
using any prior human knowledge. This is the start of our research and we will
expand to more complex problems in near future.

##  2 Related Work

The idea of automating scientific research activities dates back to the early
days of computer science. An article on _Science_ in 2009 [8] provides a great
overview on the early explorations. Also in 2009 a “Robot Scientist” named
Adam was released [9]. The authors developed specialized hardware for
conducting basic experiments, such as tracking yeast growth with varying gene
deletions and metabolites. This was paired with logic programming software for
selecting experiments. The software keeps track of various hypotheses and
chooses experiments likely to refute many of them at once. These experiments
are automatically performed, and their results guide the next experiment’s
selection. Adam effectively identified the functions of multiple genes,
requiring fewer experiments compared to other experiment-selection methods
like cost-based choices. [10] presents a research that utilizes special
hardwares to automatically learn the effects of different drugs upon the
distribution of different proteins within mammalian cells.

Very recently a breakthrough was brought by DeepMind [6], in which the authors
created a large language model that learned geometry on one billion generated
problems, in order to discover geometry properties, and train itself to prove
these properties. The model was tested on 30 IMO geometry and got 25 of them
correct, which outperforms the majority of IMO participants. This is the first
time a neutral network model learns to master a discipline of science on its
own, and it will not be surprising if the same methodology can be extended to
other disciplines such as number theory and combinatorics.

Our goal is to let AI make scientific discoveries on its own, and we start
from a simple but important problem raised hundreds of years ago: Integrals.
We are following a similar methodology as [6] to let a model learn how to
compute integrals of functions on its own.

On the other hand, there have been studies on training deep learning models to
infer the integrals of functions. In [11] the authors treated the problem of
inferring the integral of a function as a translation problem, and trained a
traditional transformer[12] using pairs of randomly generated functions and
their integrals. This is very different from our study, as it was intended to
teach a model how to infer the integral of a function by providing many
examples, while we are trying to let models discover how to compute integrals
on its own.

##  3 Experiment Design and Model Training

As mentioned in Section 1, our experiment is consisted of four steps:

  * •

Collect a large set of functions (including randomly generated ones), and for
each function $f$, compute the integral defined as area under the curve (i.e.,
$\int_{0}^{x}f(t)dt$).

  * •

For each function $f$, use a symbolic regression model to infer the function
for its integral: $g(x)=\int_{0}^{x}f(t)dt$. For example, if $f(x)=x^{2}$,
then ideally the symbolic regression model should infer that
$g(x)=\frac{x^{3}}{3}$.

  * •

Train a language model which directly infers the mathematical expression of
$g(x)$ from $f(x)$.

  * •

Try to identify simple rules of integral for each type of functions, such as
$\int_{0}^{x}t^{n}dt=\frac{x^{n+1}}{n+1}$.

We will describe each step in details in the following subsections.

###  3.1 Dataset of Functions

In this study we consider functions that are polynomial, trigonometric, or
exponential.

For polynomial functions, we get 10K of them from the _parametric functions_
data in AMPS dataset[13] (under folder _mathematica/algebra_) that can be
parsed by SymPy[14]. Then we randomize each function by randomly increase,
decrease, multiply or divide each constant (including coefficient and
exponent) by an integer between 1 and 3, without changing its sign. By
combining the original functions and the randomized functions, we get a set of
58K polynomial functions.

Then we randomly generate trigonometric or exponential functions such as
$a\cdot\sin(bx)+c$ or $a\cdot e^{bx}+c$, where $a$, $b$ and $c$ are random
constants between -10 and 10. 114K functions are generated in this way.

Obviously we could have included more types of functions. But we feel the
above set is sufficient for studying how deep learning can discover how to
calculate integrals.

For each function $f$, we compute the integral defined as area under the curve
$\int_{0}^{x}f(t)dt$. This is done using the Trapezoidal rule (by
numpy.trapz), which is simply a cumulative sum of a function over a certain
range on the $x$-axis.

###  3.2 Symbolic Regression

Symbolic regression is the procedure of inferring the mathematical expression
of a function given a sample of its values. Traditionally the standard
approach are evolutionary programs such as genetic programming
[15][16][17][18], despite their high computational cost and limited
performance. In recent years it has been found that deep neural networks
outperforms the traditional approaches [7][19][20][21].

Here we take the approach in [7] as an example. It converts a symbolic
regression problem as a sequence-to-sequence learning problem. The input to
the problem is a sample of points of the target function:
$(\vec{x_{1}},y_{1}),...,(\vec{x_{n}},y_{n})$. Each dimension of each input
point is converted into a few tokens, representing the index of the dimension
and the value. The output, which is the mathematical expression of the
function, is encoded as a sequence of tokens, each being a variable, operator,
function, or value. Then the symbolic regression problem is converted into a
sequence-to-sequence prediction problem, which tries to infer the output
sequence from the input sequence. The authors adopted a classic transformer
architecture [12], with an encoder and an decoder, each having 8-layers.
Because the transformer is a language model and is not capable of computing
function values, the constants in the output mathematical expression can
deviate from the actual values. Therefore, the last step of the approach is to
use gradient descent to optimize the constants in the output expression, until
it reaches a local minimum. There is a certain degree of randomization in the
process of inferring the output expression, i.e., generating each token of the
output sequence, and the degree of randomness can be controlled by the
_temperature_ parameter.

In order to train the model for predicting the output sequence given the input
sequence, hundreds of thousands of randomly generated functions were used in
the training process. The authors reported reasonable high accuracy ($>80\%$)
for symbolic regression problems with mathematical expressions of limited
complexity.

After comparing various approaches presented in [22] and that in [7], we found
that the approach in [7] has non-inferior accuracy but is much more time
efficient than the others, with a stable latency of about one second. This is
probably because it mainly involves a sequence-to-sequence inference using a
transformer model. In contrast, the latency of all other approaches we tested
varies greatly from case to case, from a split second to minutes. Since we
need to run symbolic regression on hundreds of thousands of cases, a low and
stable latency is key to our study, and thus we choose the approach in [7].

###  3.3 Model Training

After inferring the mathematically expression of the integral of each function
$f$ (defined as area under the curve $g(x)=\int_{0}^{x}f(t)dt$), we train a
model that directly infers the mathematical expresion of $g$ from that of $f$.
For example, given input $4x^{2}+2$, the model should output $1.33x^{3}+2x$.

Please note that we use floating point number $1.33$ instead of $\frac{4}{3}$
for the following reason. A symbolic regression model usually outputs floating
point numbers instead of fractions as coefficients, because it first
determines the structure of the mathematical expression with approximate
coefficients, and then uses gradient descent or genetic algorithm to find
optimal values for the coefficients. Therefore we have to use floating point
numbers instead of fractions.

When training a language model to infer the integral of a function, we need to
convert the mathematical expressions into sequences of tokens, and thus we
need to limit the number of digits to prevent a number from being converted
into too many tokens. For example, if there is a coefficient
“0.1428571428571”, it will be tokenized into the following tokens by GPT-Neo
[23]: “0”, “.”, “14”, “28”, “57”, “14”, “28”, “571”. In this study we choose
to retain two decimal places, which avoids the complexity introduced by such
lengthy representation of numbers. Our experiments show that keeping two
decimal places does not hurt its capability to infer the basic rules of
integrals.

We test two most popular open-source model architectures. The first one is
GPT-Neo[23]111https://github.com/EleutherAI/gpt-neo , which is a most popular
open-source implementation of the GPT-3 model (with much smaller model size).
The second is Google’s Flan T5 [24], a most popular open-source sequence-to-
sequence model.

Figure 1: An illustration of how an autoregressive model like GPT works. It
predicts each token based on all previous tokens (including a Begin-Of-
Sentence token at the beginning).

GPT-Neo is an autoregressive model, which simply predicts the next token based
on all previous tokens. Figure 1 illustrates how it works. In the transformer
architecture, each input token is initially converted into an embedding,
representing the token in a numerical form. This token embedding is then
combined with a positional embedding, which encodes the position of the token
in the sequence. The resultant combined embedding serves as the input to the
first layer of the transformer.

Within each layer of the transformer, a core component is the Multihead Self-
Attention (MSA) mechanism, which operates on the embeddings (or hidden states
from previous layers) by computing interactions between each pair of positions
in the sequence. Essentially, it allows each position to consider information
from all other positions, facilitating a comprehensive understanding of the
entire sequence.

However, to prevent the model from using future information during training (a
process known as "looking ahead"), a triangular mask is applied in the self-
attention mechanism. This mask ensures that the computation for a given
position can only utilize information from preceding positions, preserving the
sequential nature of the data.

Finally, at the output stage of the transformer, it’s the hidden state
corresponding to the last position that is used for predicting the next token
in the sequence. This approach allows the model to generate predictions based
on the entire context provided by the input sequence up to that point.

Figure 2: An illustration of a sequence-to-sequence prediction, which
translates an English sentence into French.

Flan T5 stands out as a sequence-to-sequence model, as illustrated in Figure
2. This model operates by first processing a sequence of tokens, transforming
them into embeddings which are then input into an encoder. This encoder is
essentially a stack of transformer layers, typically ranging from 12 to 48 in
number. In the encoder’s Multihead Self-Attention component within each layer,
there is no triangular mask. Consequently, each input token has the potential
to impact the hidden state at any position across the sequence.

Parallel to the encoder, the decoder in Flan T5 is also based on the
transformer architecture, generally mirroring the encoder in terms of the
number of layers. In each layer of the decoder, the Multihead Self-Attention
mechanism functions to calculate interactions between the hidden states of the
encoder and those of the decoder. This process is critical for aligning and
translating the sequence from the input to the output.

The decoder operates in an autoregressive manner, similar to the mechanism
depicted in Figure 1. It sequentially predicts each next token, relying on the
collective context provided by all previously predicted tokens. This approach
allows a sequence-to-sequence model to generate output sequences that are
coherent based on the input sequence processed by the encoder.

We tried both GPT-Neo and Flan-T5 to train models to predict the integral of a
function. Given a dataset of more than 100K functions as described in Section
3.1, for each function $f$, we compute its integral (defined as area under the
curve) $g(x)=\int_{0}^{x}f(t)dt$ using the Trapezoidal rule. Please note that
so far we only know the function values of $g(x)$ over a range on the
$x$-axis, not the mathematical expression. Then we use a symbolic regression
model to fit the above function values, which leads to a function
$g^{\prime}(x)$. The symbolic regression model is not always accurate, and we
ignore a case if

$\frac{\int_{0}^{T}|g(x)-g^{\prime}(x)|dx}{T\cdot
max(max(|g(x)|),max(|g^{\prime}(x)|))}>\epsilon$ (1)

In equation 1, $[0,T]$ is the range of our integral, and we calculate average
relative difference between $g(x)$ and $g^{\prime}(x)$, which is the average
absolute value of the gap between $g(x)$ and $g^{\prime}(x)$ divided by the
max absolute value of either function. If this average relative difference is
larger than our threshold $\epsilon$, we ignore this case because the symbolic
regression has probably been incorrect. Since we only have the values of
$g(x)$ on a large sample of points, as computed using the Trapezoidal rule, we
simply compute the average of $|g(x)-g^{\prime}(x)|$ over these points.

For all the remaining functions with integral functions $g^{\prime}(x)$
inferred, we consider the original $f$ and the inferred $g^{\prime}(x)$ as a
training example. All the training examples are split into training,
validation, and testing set. A model is trained on the training set, and we
tune the hyper parameters of the model, such as learning rate and number of
epochs trained, in order to minimize the loss on the validation set. Then we
test the model on the testing set and report the performance.

###  3.4 Discovering Rules of Integral

Every one with knowledge of calculus knows the basic rules of integral, such
as $\int x^{n}dx=\frac{1}{n+1}x^{n+1}+C$, and $\int sin(x)dx=-cos(x)+C$. In
this study we try to design a way so that these rules can be discovered
automatically. Consider the deep learning models described in Section 3.3.
They can infer the integral function from the original function (e.g., output
$0.2x^{5}$ given input $f(x)=x^{4}$), but cannot output the more generic rules
of integrals in a human-understandable manner. In this section we will discuss
how to produce these basic rules explicitly.

Here we only aim at finding very simple rules. Take polynomials as an example.
If $f(x)$ is a polynomial function, each coefficient of the integral function
$g(x)$ can be inferred from one or a few coefficients of $f(x)$ combined with
basic arithmetic operations. For example, if we consider quadratic functions,
$f(x)=a_{2}x^{2}+a_{1}x+a_{0}$, and its integral
$g(x)=a_{3}^{\prime}x^{3}+a_{2}^{\prime}x^{2}+a_{1}^{\prime}x+a_{0}$, the
basic rules would be $a_{3}^{\prime}=\frac{a_{2}}{3}$,
$a_{2}^{\prime}=\frac{a_{1}}{2}x^{2}$, $a_{1}^{\prime}=a_{0}$, and
$a_{0}^{\prime}=0$222Since we define integrals as area under the curve
starting from $x=0$, $g(0)$ is always zero.. We can also go one step further
to produce the rule for all polynomials: $a_{n}^{\prime}=\frac{1}{n}a_{n-1}$.

We design a _rule search system_ that automatically searches for such simple
rules that determines each coefficient in the integral function $g(x)$. Given
a large set of functions of the same form (e.g., $a_{4}sin(a_{3}x)+a_{0}$) It
first selects the most common form of their integral functions
($a_{4}^{\prime}cos(a_{3}^{\prime}x)+a_{1}^{\prime}x+a_{0}^{\prime}$), and
considers that to the the form of integral function $g(x)$. For each target
coefficient (i.e., coefficient in $g(x)$), it first runs a linear regression
between the target coefficient $a_{k}^{\prime}$ and each coefficient $a_{i}$
in $f(x)$ (also with $\frac{1}{a_{i}}$). If any of the linear regressions has
$|R^{2}|>0.95$, it considers that to be the formula for this target
coefficient. If no such linear regression is found, the _rule search system_
tries every pair of coefficients $a_{i}$ and $a_{j}$, combined with every
possible arithmetic operator (such as $a_{i}\cdot a_{j}$,
$\frac{a_{i}}{a_{j}}$, or $a_{i}+c\cdot a_{j}$ where $c$ is any constant). If
a linear regression between the target coefficient and the pair of
coefficients combined has $|R^{2}|>0.95$, it stops searching and declares
success. Otherwise it keeps searching for triples of coefficients, and stops
there (either succeeds or fails).

Our experiments show that the above system successfully discovers formula for
integrating functions in our dataset. For example, it discovers that for
$f(x)=a_{4}exp(a_{3}x)+a_{0}$, its integral is
$\int_{0}^{x}f(t)dt=\frac{a_{4}}{a_{3}}exp(a_{3}x)+a_{0}x-\frac{a_{4}}{a_{3}}$.
Similarly it discovers that for $f(x)=a_{4}sin(a_{3}x)+a_{0}$, its integral is
$\int_{0}^{x}f(t)dt=-\frac{a_{4}}{a_{3}}cos(a_{3}x)+a_{0}x+\frac{a_{4}}{a_{3}}$.

For polynomials, it discovers

$\begin{split}\int_{0}^{x}a_{0}+a_{1}t+a_{2}t^{2}+a_{3}t^{3}+a_{4}t^{4}+a_{5}t^{5}+a_{6}t^{6}dt=\\\
a_{0}t+0.5a_{1}t^{2}+0.333a_{2}t^{3}+0.25a_{3}t^{4}+0.2a_{4}t^{5}+0.167a_{5}t^{6}+0.143a_{6}t^{7}\end{split}$

It also discovers that the new coefficient for $k$-th power
$a_{k}^{\prime}=\frac{a_{k-1}}{k}$.

##  4 Experiment Results

We use the GPT-Neo model[23] of 350M trainable parameters, downloaded from
Huggingface333https://huggingface.co/xhyi/PT_GPTNEO350_ATG, and the
Flan-T5-large model [24], which has 783M parameters, also downloaded from
Huggingface444https://huggingface.co/google/flan-t5-base. All experiments are
done on a machine with an A6000 GPU, with Ubuntu 18.04, CUDA 11.7 and Pytorch
2.0.0. Float32 is used because numerical precision is key to our approach.

###  4.1 Function Dataset

As described in Section 3.1, we get 9,024 polynomial functions from the
_parametric functions_ dataset in AMPS dataset[13] that can be parsed by
SymPy. Then we randomize each of them by randomly increase, decrease, multiply
or divide each constant (including coefficient and exponent) by an integer
between 1 and 3 without changing its sign. In this way we get 58,512
polynomial functions in total.

Then we randomly generate trigonometric or exponential functions such as
$a\cdot\sin(bx)+c$, $a\cdot\cos(bx)+c$ and $a\cdot e^{bx}+c$, where $a$, $b$
and $c$ are random constants between -10 and 10. 113,825 functions were
generated in this way.

Obviously we could have included more types of functions. But we feel the
above set is sufficient for studying how deep learning can discover the way to
calculate integrals.

###  4.2 Symbolic Regression

For each function $f$ in our function dataset, we compute the integral defined
as area under the curve $g(x)=\int_{0}^{x}f(t)dt$. We use numerical methods to
compute the value of $g(x)$ on the range $[0,T]$, using the Trapezoidal
rule555Using function numpy.trapz. which is simply a cumulative sum of a
function over $[0,T]$ on the $x$-axis. We set $T=4$, and compute values on
10,000 points for each function.

This gives us 10,000 pairs of $(x,g(x))$ for each function. We apply the
symbolic regression model [7] mentioned in Section 3.2, which takes a number
of $(x,y)$ pairs as input and outputs a function that is inferred by the
underlying deep learning model, whose coefficients are adjusted by gradient
descent to minimize the error against the input points. We denote the output
function of symbolic regression as $g^{\prime}(x)$.

The symbolic regression model is not always accurate, and we need to filter
out cases where $g^{\prime}(x)$ is obviously incorrect. As discussed in the
end of Section 3.3, we compute the average relative absolute difference
between $g(x)$ and $g^{\prime}(x)$ over the 10,000 points, and ignore a case
if the difference is greater than $\epsilon=0.01$. After this filtering we
have 49,109 polynomial functions left (out of 58,512), and 63,409 non-
polynomial functions left (out of 113,825). Out of these 112.5K examples, We
use 90K as our training set, and 22.5K as our evaluation set which is used to
tune the hyper-parameters. We generated additional 10K polynomial and 10K non-
polynomial functions in the same way as our testing dataset.

When running symbolic regression, it is possible that given the data points of
$cos(x)$, a symbolic regression model outputs $sin(\frac{\pi}{2}-x)$. It is
also possible that given the data points of $a\cdot e^{bx}$, the model outputs
$a^{\prime}\cdot e^{bx-c}$, where $a^{\prime}=a\cdot e^{c}$. Fortunately the
symbolic regression approach we choose prefers simpler outputs, and thus in
the majority of cases we do not observe the above issue.

Please note that a $g^{\prime}(x)$ can still be incorrect even if its relative
error is smaller than $\epsilon$. Since our goal is to use deep learning to
learn to infer the integral of functions purely on its own (without human
supervision), we have to live with the errors of symbolic regression. Later we
will show that our models are more accurate than symbolic regression, and are
much faster because only it only needs to infer a sequence of tokens.

###  4.3 Training Models for Integral

As described above, we have acquired 112.5K examples, each having a function
$f(x)$ and its integral $g(x)=\int_{0}^{x}f(t)dt$. For example,
$f(x)=2x^{2}+3.55x-1.4$, and $g(x)=0.67x^{3}+1.78x^{2}-1.4x$ is the output of
symbolic regression. Please note that the raw output of symbolic regression
can have many digits (e.g., $1.77501292$) and we round all coefficients to two
decimal places in our study.

Section 3.3 described the two models we use in our study: An autoregressive
model called GPT-Neo, and a sequence-to-sequence model called Flan-T5. We
train both models to learn to predict the integral $g(x)$ by taking the
original function $f(x)$ as input. For GPT-Neo, we need to convert each pair
of $f(x)$ and $g(x)$ into a string, which will be split into tokens and fed
into the GPT-Neo model to learn to predict the next token based on previous
tokens. We use SymPy [14] to convert a math expression into a string. For the
above example, it is converted into

“$2*x**2+3.55*x-1.4\>entail\>0.67*x**3+1.78*x**2-1.4*x$”

Please note we use a special token "entail" to separate the original function
and its integral. At inference time, we can feed “$2*x**2+3.55*x-1.4\>entail$”
into the model, and it should output the integral, one token at a time.

Unlike GPT-Neo, Flan-T5 is a sequence-to-sequence model, and we simply use the
original function (such as “$2*x**2+3.55*x-1.4$”) as the input sequence, and
the integral as the output sequence (such as “$0.67*x**3+1.78*x**2-1.4*x$”).

####  4.3.1 Accuracy Evaluation

Our models may not always produce the correct expression for integral. It
might be due to the model’s capability, or due to errors in its training data
generated by symbolic regression. Given a function $f$ and its integral
$g(x)=\int_{0}^{x}f(t)dt$ (rounded to two decimal places), suppose a model’s
output is $\hat{g}(x)$. We say $\hat{g}(x)$ is correct if $\hat{g}(x)=g(x)$.

Sometimes $\hat{g}(x)$ only differs a little from $g(x)$, which is probably
due to rounding errors. Consider our example mentioned above:
$f(x)=2x^{2}+3.55x-1.4$ and $g(x)=0.67x^{3}+1.78x^{2}-1.4x$. Our model may
produce $\hat{g}(x)=0.67x^{3}+1.77x^{2}-1.4x$, since symbolic regression may
produce such a result which is almost equally good as $g(x)$. We define the
error between $g(x)$ and $\hat{g}(x)$ in the same way as we define the error
of symbolic regression:

$RelativeDiff(g(x),\hat{g}(x))=\frac{\int_{0}^{T}|g(x)-\hat{g}(x)|dx}{T\cdot
max(max(|g(x)|),max(|\hat{g}(x)|))}$ (2)

We say $\hat{g}(x)$ is approximately correct if
$RelativeDiff(g(x),\hat{g}(x))<\epsilon$, where $\epsilon=0.01$. In the above
example, $\hat{g}(x)=0.67x^{3}+1.77x^{2}-1.4x$ is not correct but is
approximately correct.

####  4.3.2 Pre-training

Although there are open-source pre-trained version of each model available on
the web, we do not use them because the pre-training may involve knowledge
about calculus. Instead, we pre-train each model using the _algebra_ and
_number theory_ datasets in the AMPS$\rightarrow$Mathematica dataset [13],
which contain no knowledge of calculus. These datasets contain 1.99M math
problems and answers, mostly at grade-school level. The pre-training should
teach each model about basic math knowledge on algebra and number theory.

Table 1 shows how GPT-Neo’s accuracy improves with pre-training. We can see
that pre-training improves the percentage of correct from 74.1% to 80.2%, and
that of approximately correct from 85.9% to 95.3%. The accuracy is about 2%
lower than that with the open-source model, probably because the latter has
been trained with some knowledge about calculus.

Table 1: GPT-Neo’s Performances with various amounts of pre-training. In each case the model is first trained with pre-trainined data, and then trained five epochs on the our dataset for inferring integral of functions Pre-Training | %-Correct | %-Approx.-Correct
---|---|---
Randomly Initialized | 74.1% | 85.9%
After 2 epochs | 78.4% | 93.0%
After 5 epochs | 80.2% | 95.2%
After 8 epochs | 80.2% | 95.3%
From HuggingFace111Downloaded from https://huggingface.co/xhyi/PT_GPTNEO350_AT | 82.4% | 97.3%
|  |

####  4.3.3 Model Accuracy

In this section we report the accuracy of our models in inferring integral of
a function. As mentioned in Section 4.2, we have a training set of 90K
examples, a validation set of 22.5K, and a testing set of 20K (half
polynomials and half non-polynomials). We train our models using the training
set, based on the pre-trained model. We use the validation set to tune the
hyper-parameters, and we report accuracy on the testing set.

For Flan-T5, we found the validation loss stabilized in the fifth epoch, and
thus we trained the model for five epochs. We used learning rate of $3\times
10^{-4}$ and batch size of 8. Figure 3 shows its accuracy on inferring
integrals of polynomial functions, and how it changes w.r.t. number of
examples trained (90K $\times$ 5 epochs = 450K examples in total). We can see
that the %-Correct grows from 46.2% to 67.7%, and %-Approx.-Correct grows from
64.0% to 85.7%.

Figure 4 shows Flan-T5’s accuracy on inferring integrals of non-polynomial
functions. The %-Correct grows from 45.1% to 64.6%, and %-Approx.-Correct
grows from 52.2% to 75.7%.

We can see Flan-T5 is not highly accurate, and we find that GPT-Neo has much
higher accuracy, as shown below.

$0$$100$$200$$300$$400$$0$$20$$40$$60$$80$$100$Examples Trained
(Thousands)Accuracy (%)%-Correct%-Approx.-Correct

Figure 3: The %-Correct and %-Approximately-Correct of Flan-T5 in inferring
integrals of polynomial functions.

$0$$100$$200$$300$$400$$0$$20$$40$$60$$80$$100$Examples Trained
(Thousands)Accuracy (%)%-Correct%-Approx.-Correct

Figure 4: The %-Correct and %-Approximately-Correct of Flan-T5 in inferring
integrals of non-polynomial functions.

For GPT-Neo, we also trained the model for five epochs, using learning rate of
$2\times 10^{-5}$ and batch size of 16. For polynomials, %-Correct improves
from 57.0% to 73.3% over the course of training, and %-Approximately-Correct
(i.e., no more than 1% average relative error, which is probably due to
rounding errors in training data caused by symbolic regression) improves from
77.6% to 95.7%.

For non-polynomials, %-Correct improves from 67.4% to 87.1% over the course of
training, and %-Approximately-Correct improves from 74.9% to 94.8%.

$0$$100$$200$$300$$400$$0$$20$$40$$60$$80$$100$Examples Trained
(Thousands)Accuracy (%)%-Correct%-Approx.-Correct

Figure 5: The %-Correct and %-Approximately-Correct of GPT-Neo in inferring
integrals of polynomial functions.

$0$$100$$200$$300$$400$$0$$20$$40$$60$$80$$100$Examples Trained
(Thousands)Accuracy (%)%-Correct%-Approx.-Correct

Figure 6: The %-Correct and %-Approximately-Correct of GPT-Neo in inferring integrals of non-polynomial functions. Table 2: Comparing Performances of GPT-Neo, Flan-T5, and Symbolic Regression in computing integrals | GPT-Neo | Flan-T5 | Symbolic Regression
---|---|---|---
%-Correct Polynomials | 73.3% | 67.7% | 65.7%
%-Approx-Correct Polynomials | 95.7% | 85.7% | 90.6%
%-Correct Non-polynomials | 87.1% | 64.6% | 46.4%
%-Approx-Correct Non-polynomials | 94.8% | 75.7% | 71.1%
Avg. Time (sec) | 0.722 | 0.151 | 3.17
|  |  |

Table 2 compares the performances of GPT-Neo, Flan-T5, and the Symbolic
Regression model in [7]. Please note that GPT-Neo and Flan-T5 directly infer
the integral function from the original function, without any numeric
computation. In contrast, in order to use the Symbolic Regression model to
infer the integral of a function, we need to first use the Trapezoidal rule to
compute the values of the integral function, and then use the model to infer
the math expression of the function.

From Table 2 we can see that GPT-Neo is the most accurate among the three. As
mentioned before, the training data of GPT-Neo and Flan-T5 contain some
rounding errors, because the training data was generated using symbolic
regression which is not 100% precise. Such small errors in the training data
leads to some small errors (usually $<$1%) in the output integral functions of
GPT-Neo and Flan-T5. We can see that although GPT-Neo produces exactly the
expected integral function in about 80% of cases, it has 95%+ accuracy if we
allow such small errors ($<$1%).

###  4.4 Finding Rules for Integrals

####  4.4.1 Rules of integral for polynomials

As described in Section 3.4, we aim at finding very simple rules for inferring
of integrals of each type of functions. For example, if we consider quadratic
functions as a type, the integral of function $f(x)=a_{2}x^{2}+a_{1}x+a_{0}$
is $g(x)=a_{3}^{\prime}x^{3}+a_{2}^{\prime}x^{2}+a_{1}^{\prime}x+a_{0}$, the
basic rules would be $a_{3}^{\prime}=\frac{a_{2}}{3}$,
$a_{2}^{\prime}=\frac{a_{1}}{2}x^{2}$, $a_{1}^{\prime}=a_{0}$, and
$a_{0}^{\prime}=0$666$a_{0}^{\prime}=0$ because we define integrals as area
under the curve starting from $x=0$..

Please note that we compute the integral of a function by first applying
Trapezoidal rule to compute the integral function’s values, and then using the
symbolic regression model in [7] to find the math expression, and excluding
all cases where the average relative error is greater than 1%. Therefore, the
discovery of rules of integrals does not depend on our ability to train models
to infer integral functions.

In this study, we design a _rule search system_ , which tries every possible
combination of single, double (and more) coefficients, combined together with
arithmetic operators, and checks if that can be used to infer a coefficient in
the integral function $g(x)$. In order to find a rule to infer a particular
coefficient $a_{k}^{\prime}$ in $g(x)$, the _rule search system_ starts from
trying each individual coefficient $a_{i}$ (and also $\frac{1}{a_{i}}$) in
$f(x)$ by running linear regression with $a_{i}$ and $a_{k}^{\prime}$. If it
finds an $a_{i}$ that leads to a high $R^{2}$ in the linear regression
($|R^{2}|>0.95$), it stops and declares success. Otherwise it tries each pair
of coefficients $a_{i}$ and $a_{j}$ in $f(x)$, by combining them using each
arithmetic operation. If no success with pairs of coefficients, it tries
triples, and stops there.

Table 3: $R^{2}$ of linear regression between individual coefficients of
$f(x)$ ($a_{0},...,a_{6}$) and those of $g(x)$
($a_{1}^{\prime},...,a_{7}^{\prime}$), for polynomial functions

$R^{2}$ | $a_{0}$ | $a_{1}$ | $a_{2}$ | $a_{3}$ | $a_{4}$ | $a_{5}$ | $a_{6}$
---|---|---|---|---|---|---|---
$a_{1}^{\prime}$ | 1.000 | 0.006 | 0.177 | 0.036 | 0.095 | -0.021 | -0.031
$a_{2}^{\prime}$ | 0.006 | 1.000 | 0.002 | 0.137 | 0.010 | 0.005 | -0.005
$a_{3}^{\prime}$ | 0.176 | -0.001 | 1.000 | -0.045 | -0.115 | -0.027 | -0.043
$a_{4}^{\prime}$ | 0.036 | 0.142 | -0.033 | 1.000 | -0.016 | -0.006 | -0.010
$a_{5}^{\prime}$ | 0.095 | 0.009 | -0.114 | -0.023 | 1.000 | -0.015 | -0.023
$a_{6}^{\prime}$ | -0.020 | 0.005 | -0.027 | -0.006 | -0.015 | 1.000 | -0.020
$a_{7}^{\prime}$ | -0.031 | -0.002 | 0.039 | -0.006 | -0.022 | 0.020 | 1.000
|  |  |  |  |  |  |

66footnotetext: Note: $a_{0}^{\prime}\approx 0$ in almost all cases and thus
does not require further studies.

Table 3 shows the $R^{2}$ of linear regression between individual coefficients
of $f(x)$ and those of $g(x)=\int_{0}^{x}f(t)dt$, where $f(x)$ is a polynomial
function and $a_{i}$ is the coefficient for $x^{i}$. We can see that each
coefficient of $g(x)$ ($a_{1}^{\prime},...,a_{7}^{\prime}$) has $R^{2}$ of 1.0
with a particular coefficient of $f(x)$ ($a_{0},...,a_{6}$). To be more
specific, $a_{k}^{\prime}$ has perfect correlation with $a_{k-1}$, and has
little correlation with other coefficients of $f(x)$. This shows that each
$a_{k}^{\prime}$ can be perfectly predicted by $a_{k-1}$.

Table 4: Slope and intercept of linear regression between each coefficient $a_{k}^{\prime}$ of $g(x)$ and the corresponding coefficient in $f(x)$ that is correlated with $a_{k}^{\prime}$ (i.e. $a_{k-1}$), for polynomial functions  | Slope | Intercept |  |  |  |
---|---|---|---|---|---|---
$a_{1}^{\prime}$ | 1.000 | -0.003 |  |  |  |
$a_{2}^{\prime}$ | 0.500 | -0.002 |  |  |  |
$a_{3}^{\prime}$ | 0.333 | -0.003 |  |  |  |
$a_{4}^{\prime}$ | 0.250 | 0.000 |  |  |  |
$a_{5}^{\prime}$ | 0.200 | 0.000 |  |  |  |
$a_{6}^{\prime}$ | 0.167 | 0.000 |  |  |  |
$a_{7}^{\prime}$ | 0.143 | 0.000 |  |  |  |
|  |  |  |  |  |

Table 4 shows the slope and intercept of the linear regression between each
coefficient $a_{k}^{\prime}$ of $g(x)$ and the corresponding coefficient in
$f(x)$ that is correlated with $a_{k}^{\prime}$ (i.e., $|R^{2}|$ is close to
1.0). We can see that the intercept is close to zero, and therefore the
formula for integrating polynomials (up to 6th power) can be written as

$\begin{split}\int_{0}^{x}a_{0}+a_{1}t+a_{2}t^{2}+a_{3}t^{3}+a_{4}t^{4}+a_{5}t^{5}+a_{6}t^{6}dt=\\\
a_{0}t+0.5a_{1}t^{2}+0.333a_{2}t^{3}+0.25a_{3}t^{4}+0.2a_{4}t^{5}+0.167a_{5}t^{6}+0.143a_{6}t^{7}\end{split}$
(3)

Now we let our rule search system find the connection between the slope
$\beta$ and the power $p$ (e.g., $a_{i}^{\prime}$ is associated with power
$i$), because they are the only two meaningful variables left. It first runs
linear regression between $p$ and $\beta$, and gets $R^{2}$ of 0.735, which is
not high enough to confirm their connection. Then it tries to find the
relationship between $\beta$ and a simple transform of $p$ using any
arithmetic operator, and finds that $\frac{1}{p}$ has perfect correlation
($R^{2}=1.00$) with $\beta$, and $\beta=1.00\cdot\frac{1}{p}+0.001$. Our rule
search system discards any constants smaller than 0.01 because they are likely
caused by rounding errors. Thus it concludes that $\beta=\frac{1}{p}$.

Combining with the results in 3, it concludes that the coefficients of $g(x)$
can be inferred as

$a_{k}^{\prime}=\frac{1}{k}a_{k-1},$ (4)

which means that

$\int_{0}^{x}a_{k-1}t^{k-1}dt=\frac{a_{k-1}}{k}x^{k}$ (5)

####  4.4.2 Rules of integral for exp function

Now let us look at non-polynomial functions. Since we studied three types of
non-polynomial functions ($exp$, $sin$ and $cos$), we need to study their
integral rules separately.

We start by applying our _rule search system_ on the 15.3K cases involving
$exp$ functions, which take the form of $a_{4}exp(a_{3}x+a_{2})+a_{1}x+a_{0}$,
where $a_{2}$ and $a_{1}$ are both zero in our generated functions. The _rule
search system_ found that among their integral functions, 99.9% of them take
the form of
$a_{4}^{\prime}exp(a_{3}^{\prime}x+a_{2}^{\prime})+a_{1}^{\prime}x+a_{0}^{\prime}$,
and concluded that this is the form of $exp$ function’s integral. It was also
found that $a_{2}^{\prime}$ are mostly zero and thus require no further
investigation. The correlation (indicated by $R^{2}$ of linear regression) of
the remaining coefficients is shown in Table 5.

Table 5: $R^{2}$ of linear regression between individual coefficients of
$f(x)$ ($a_{0},...,a_{4}$) and those of the integral function $g(x)$
($a_{0}^{\prime},...,a_{4}^{\prime}$) , where
$f(x)=a_{4}exp(a_{3}x+a_{2})+a_{1}x+a_{0}$

$R^{2}$ | $a_{0}$ | $a_{1}$ | $a_{2}$ | $a_{3}$ | $a_{4}$
---|---|---|---|---|---
$a_{0}^{\prime}$ | 0.121 | N/A | N/A111$a_{1}$ and $a_{2}$ are all zeros and linear regression could not be performed. | -0.067 | 0.707
$a_{1}^{\prime}$ | 1.000 | N/A | N/A | -0.027 | 0.133
$a_{3}^{\prime}$ | -0.027 | N/A | N/A | 1.000 | -0.013
$a_{4}^{\prime}$ | -0.121 | N/A | N/A | 0.069 | -0.706
|  |  |  |  |

66footnotetext: Note: $a_{2}^{\prime}=0$ in most cases and thus does not
require further studies.

From Table 5 one can see that $a_{1}^{\prime}$ is perfectly correlated with
$a_{0}$ and $a_{3}^{\prime}$ perfectly correlated with $a_{3}$, which are
consistent with our expectation. Based on the results of linear regressions,
it is found that $a_{1}^{\prime}\approx a_{0}$ (slope 1.000 and intercept
0.000), and $a_{3}^{\prime}\approx a_{3}$ (slope 1.000 and intercept 0.000).
But there is no conclusion drawn for $a_{0}^{\prime}$ and $a_{4}^{\prime}$,
and the _rule search system_ investigates pairs of coefficients among
$a_{0},...,a_{4}$, combined with an arithmetic operation.

Table 6: The maximum $R^{2}$ of linear regression between $a_{4}^{\prime}$ (of integral of a function involving $exp$) and any pair of coefficients of the original function, combined by an arithmetic operator. We also show the slope and intercept of linear regression.  Operator |  $R^{2}$ with max abs. value | Slope | Intercept | First Coefficient | Second Coefficient
---|---|---|---|---|---
Addition111A multi-variate linear regression is performed to use the pairs of coefficients to predict $a_{4}^{\prime}$, which allows the pair of coefficients to be weighted summed. | -0.706222A linear regression using $a_{0},...,a_{4}$ and $a_{4}^{\prime}$ also has $|R^{2}|$ of 0.706. | -0.124 | -0.013 | $a_{4}$ | $a_{4}$
Multiplication | 0.807 | 0.080 | -0.033 | $a_{3}$ | $a_{4}$
Division | 1.000 | 1.000 | 0.001 | $a_{4}$ | $a_{3}$
|  |  |  |  |

From Table 6 it is found that the only combination of coefficients leading to
qualified result ($|R^{2}|>0.95$) is $\frac{a_{4}}{a_{3}}$. The formula found
is $a_{4}^{\prime}=1.000\frac{a_{4}}{a_{3}}+0.001$, and we can confidently
draw the conclusion of $a_{4}^{\prime}=\frac{a_{4}}{a_{3}}$777Any constant
smaller than 0.01 is discarded because it is probably caused by rounding
errors..

Table 7: The maximum $R^{2}$ of linear regression between $a_{0}^{\prime}$ (of integral of a function involving $exp$) and any pair of coefficients of the original function, combined by an arithmetic operator. We also show the slope and intercept of linear regression.  Operator |  $R^{2}$ with max abs. value | Slope | Intercept | First Coefficient | Second Coefficient
---|---|---|---|---|---
Addition111A multi-variate linear regression is performed to use the pairs of coefficients to predict $a_{4}^{\prime}$, which allows the pair of coefficients to be weighted summed. | -0.706 | -0.124 | -0.013 | $a_{4}$ | $a_{4}$
Multiplication | -0.807 | -0.080 | 0.034 | $a_{3}$ | $a_{4}$
Division | -0.999 | -0.999 | 0.000 | $a_{4}$ | $a_{3}$
|  |  |  |  |

Similarly we find out that $a_{0}^{\prime}=\frac{a_{4}}{a_{3}}$, with results
in Table 7. Combining with previous results, our _rule search system_ finds
that for $f(x)=a_{4}exp(a_{3}x)+a_{0}$, its integral is
$\int_{0}^{x}f(t)dt=\frac{a_{4}}{a_{3}}exp(a_{3}x)+a_{0}x-\frac{a_{4}}{a_{3}}$.

####  4.4.3 Rules of integral for sin function

Then we apply our _rule search system_ on the 23.6K cases involving $sin$
functions, which take the form of $a_{4}sin(a_{3}x+a_{2})+a_{1}x+a_{0}$, where
$a_{2}$ and $a_{1}$ are both zero in our generated functions. The _rule search
system_ found that among their integral functions, 93.8% of them take the form
of
$a_{4}^{\prime}cos(a_{3}^{\prime}x+a_{2}^{\prime})+a_{1}^{\prime}x+a_{0}^{\prime}$,
and concluded that this is the form of $sin$ function’s integral. It was also
found that $a_{2}^{\prime}$ are mostly zero and thus require no further
investigation. The correlation (indicated by $R^{2}$ of linear regression) of
the remaining coefficients is shown in Table 8.

Table 8: $R^{2}$ of linear regression between individual coefficients of
$f(x)$ ($a_{0},...,a_{4}$) and those of the integral function $g(x)$
($a_{0}^{\prime},...,a_{4}^{\prime}$), where $f(x)=a_{4}sin(a_{3}x)++a_{0}$

$R^{2}$ | $a_{0}$ | $a_{1}$ | $a_{2}$ | $a_{3}$ | $a_{4}$
---|---|---|---|---|---
$a_{0}^{\prime}$ | 0.070 | N/A | N/A111$a_{1}$ and $a_{2}$ are all zeros and linear regression could not be performed. | -0.036 | 0.900
$a_{1}^{\prime}$ | 1.000 | N/A | N/A | 0.056 | 0.073
$a_{3}^{\prime}$ | 0.058 | N/A | N/A | 0.995 | -0.012
$a_{4}^{\prime}$ | -0.062 | N/A | N/A | 0.032 | -0.885
|  |  |  |  |

77footnotetext: Note: $a_{2}^{\prime}=0$ in most cases and thus do not require
further studies.

From Table 8 we can see that $a_{1}^{\prime}$ is perfectly correlated with
$a_{0}$, and the same between $a_{3}^{\prime}$ and $a_{3}$, which are
consistent with our expectations. Based on the results of linear regressions,
it is found that $a_{1}^{\prime}=a_{0}$ (slope 1.000 and intercept 0.000), and
$a_{3}^{\prime}\approx a_{3}$ (slope 0.999 and intercept -0.001). Again there
is no conclusion for $a_{0}^{\prime}$ and $a_{4}^{\prime}$, which has some
correlation with $a_{4}$ but is lower than our requirement of $|R^{2}|>0.95$.
Thus the _rule search system_ needs to explore pairs of coefficients combined
with arithmetic operators.

Table 9: The maximum $R^{2}$ of linear regression between $a_{4}^{\prime}$ (of integral of $sin$ function) and any pair of coefficients of $f(x)$, combined by each operator. We also show the slope and intercept of linear regression.  Operator |  $R^{2}$ with max abs. value | Slope | Intercept | First Coefficient | Second Coefficient
---|---|---|---|---|---
Addition111A multi-variate linear regression is performed to use the pairs of coefficients to predict $a_{4}^{\prime}$, which allows the pair of coefficients to be weighted summed. | -0.894 | -0.276 | 0.037 | $a_{4}$ | $a_{3}$
Multiplication | -0.763 | -0.052 | 0.012 | $a_{3}$ | $a_{4}$
Division | 0.997 | 0.997 | -0.001 | $a_{4}$ | $a_{3}$
|  |  |  |  |

From Table 9 it is found that the only combination of coefficients leading to
qualified result ($|R^{2}|>0.95$) is $\frac{a_{4}}{a_{3}}$. The formula found
is $a_{4}^{\prime}=0.997\frac{a_{4}}{a_{3}}-0.001$. Since anything less than
0.01 is likely to be caused by rounding errors, we can draw the conclusion of
$a_{4}^{\prime}=\frac{a_{4}}{a_{3}}$.

Table 10: The maximum $R^{2}$ of linear regression between $a_{0}^{\prime}$ (of integral of $sin$ function) and any pair of coefficients of $f(x)$, combined by each operator. We also show the slope and intercept of linear regression.  Operator |  $R^{2}$ with max abs. value | Slope | Intercept | First Coefficient | Second Coefficient
---|---|---|---|---|---
Addition111A multi-variate linear regression is performed to use the pairs of coefficients to predict $a_{4}^{\prime}$, which allows the pair of coefficients to be weighted summed. | 0.900 | 0.292 | 0.016 | $a_{4}$ | $a_{3}$
Multiplication | 0.783 | 0.065 | -0.007 | $a_{3}$ | $a_{4}$
Division | 0.996 | 0.998 | -0.001 | $a_{4}$ | $a_{3}$
|  |  |  |  |

Similarly the _rule search system_ finds that
$a_{0}^{\prime}=-\frac{a_{4}}{a_{3}}$, with results in Table 10. Combining
previous results, it finds that for $f(x)=a_{4}sin(a_{3}x)+a_{0}$, its
integral is
$\int_{0}^{x}f(t)dt=-\frac{a_{4}}{a_{3}}cos(a_{3}x)+a_{0}x+\frac{a_{4}}{a_{3}}$.

Integrating functions involving $cos$ is similar to that for $sin$, and the
_rule search system_ concludes that for $f(x)=a_{4}cos(a_{3}x)+a_{0}$, its
integral is $\int_{0}^{x}f(t)dt=\frac{a_{4}}{a_{3}}sin(a_{3}x)+a_{0}x$.

##  5 Conclusions and Future Work

The recent advances in large language models have shown that AI can solve
difficult math problems or write code for programming challenges. However,
being able to solve a type of well-defined problems after being trained on a
huge corpus is very different from making discoveries in scientific research.

In this paper we explore how to use deep learning models to discover an old
but most important math tool: Integrals. We show that with a definition of
integrals (defined as area under the curve in our study), AI models can learn
to infer the integral of a function, such as inferring that
$\int_{0}^{x}t^{2}dt=\frac{x^{3}}{3}$ and
$\int_{0}^{x}ae^{bt}dt=\frac{a}{b}e^{bx}-\frac{a}{b}$. Our experiments show
that deep learning models can learn to infer the integral function using
autogressive or sequence-to-sequence models (similar to generating text or
translating from one language to another one), with high accuracy. We also
designed a _rule search system_ which can automatically discover the basic
rules of integrals such as $\int_{0}^{x}t^{n}dt=\frac{x^{n+1}}{n+1}$.

Discovering how to compute integrals is just the first step of our journey
towards automated research. In future we will explore how to let deep learning
models learn to solve more sophisticated problems on their own, and the
ultimate goal is to enable AI to assist human in making real scientific
discoveries.

## References

  *   * OpenAI [2023] OpenAI: Gpt-4 technical report. (2023) arXiv:2303.08774
  * [2] Microsoft Copilot. https://copilot.microsoft.com/ (2023)
  * Rozière et al. [2023] Rozière, B., et al.: Code llama: Open foundation models for code. (2023) arXiv:2308.12950
  * Huang et al. [2023] Huang, Y., et al.: Competition-level problems are effective llm evaluators. (2023) arXiv:2312.02143
  * Azerbayev et al. [2023] Azerbayev, Z., et al.: Llemma: An open language model for mathematics. (2023) arXiv:2310.10631
  * Trinh et al. [2024] Trinh, T.H., et al.: Solving olympiad geometry without human demonstrations. Nature 625, 476–482 (2024)
  * d’Ascoli et al. [2022] d’Ascoli, S., Kamienny, P.-A., Lample, G., Charton, F.: Deep symbolic regression for recurrent sequences. In: Proceedings of the 39th International Conference on Machine Learning, Baltimore, Maryland, USA (2022)
  * Waltz and Buchanan [2009] Waltz, D., Buchanan, B.G.: Automating science. Science 324(5923), 43–44 (2009)
  * King et al. [2009] King, R.D., et al.: The robot scientist adam. Computer 42(8), 46–54 (2009)
  * Naik et al. [2016] Naik, A.W., et al.: Active machine learning-driven experimentation to determine compound effects on protein patterns. eLife 5(e10047) (2016)
  * Lample and Charton [2019] Lample, G., Charton, F.: Deep learning for symbolic mathematics (2019) arXiv:1912.01412
  * Vaswani et al. [2017] Vaswani, A., et al.: Attention is all you need. In: Proceedings of the 31st Conference on Advances in Neural Information Processing Systems, Long Beach, California, USA (2017)
  * Hendrycks et al. [2021] Hendrycks, D., et al.: Measuring mathematical problem solving with the math dataset (2021) arXiv:2103.03874
  * Meurer et al. [2017] Meurer, A., et al.: Sympy: symbolic computing in python. PeerJ Computer Science 3, 103 (2017)
  * Augusto and Barbosa [2000] Augusto, D.A., Barbosa, H.J.: Symbolic regression via genetic programming. In: Proceedings Vol. 1. Sixth Brazilian Symposium on Neural Networks (2000)
  * Schmidt and Lipson [2009] Schmidt, M., Lipson, H.: Distilling free-form natural laws from experimental data. Science 324(5923), 81–85 (2009)
  * Murari et al. [2014] Murari, A., Peluso, E., Gelfusa, M., Lupelli, I., Lungaroni, M., Gaudio, P.: Symbolic regression via genetic programming for data driven derivation of confinement scaling laws without any assumption on their mathematical form. Plasma Physics and Controlled Fusion 57(1) (2014)
  * McKay et al. [1995] McKay, B., Willis, M.J., Barton, G.W.: Using a tree structured genetic algorithm to perform symbolic regression. In: Proceedings of First International Conference on Genetic Algorithms in Engineering Systems: Innovations and Applications (1995)
  * Sahoo et al. [2018] Sahoo, S., Lampert, C., Martius, G.: Learning equations for extrapolation and control. In: Proceedings of International Conference on Machine Learning (2018)
  * Kim et al. [2020] Kim, S., et al.: Integration of neural network-based symbolic regression in deep learning for scientific discovery. IEEE Transactions on Neural Networks and Learning Systems 32(9) (2020)
  * Petersen et al. [2019] Petersen, B.K., et al.: Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients (2019) arXiv:1912.04871
  * Landajuela et al. [2022] Landajuela, M., et al.: A unified framework for deep symbolic regression. In: Advances in Neural Information Processing Systems (2022)
  * Black et al. [2021] Black, S., Gao, L., Wang, P., Leahy, C., Biderman, S.: GPT-Neo: Large Scale Autoregressive Language Modeling with Mesh-Tensorflow. https://doi.org/10.5281/zenodo.5297715
  * Chung et al. [2022] Chung, H.W., et al.: Scaling instruction-finetuned language models (2022) arXiv:2210.11416


2022

[1,2]Jia Li

[1]State Key Laboratory of Virtual Reality Technology and Systems, School of
Computer Science and Engineering, Beihang University, Beijing, 100191, China

2]National Astronomical Observatories, Chinese Academy of Sciences, Beijing,
100101, China

# Deep Learning Eliminates Massive Dust Storms from Images of Tianwen-1

Hongyu Li<EMAIL_ADDRESS><EMAIL_ADDRESS>Xin Ren<EMAIL_ADDRESS>Long Xu<EMAIL_ADDRESS>*  [

###### Abstract

Dust storms may remarkably degrade the imaging quality of Martian orbiters and
delay the progress of mapping the global topography and geomorphology. To
address this issue, this paper presents an approach that reuses the image
dehazing knowledge obtained on Earth to resolve the dust-removal problem on
Mars. In this approach, we collect remote-sensing images captured by Tianwen-1
and manually select hundreds of clean and dusty images. Inspired by the haze
formation process on Earth, we formulate a similar visual degradation process
on clean images and synthesize dusty images sharing a similar feature
distribution with realistic dusty images. These realistic clean and synthetic
dusty image pairs are used to train a deep model that inherently encodes dust
irrelevant features and decodes them into dust-free images. Qualitative and
quantitative results show that dust storms can be effectively eliminated by
the proposed approach, leading to obviously improved topographical and
geomorphological details of Mars.

###### keywords:

Dust Storm, Tianwen-1, Mars, Image Enhancement

Almost all the fantastic stories (Robinson, 1993; Scott, 2015) of the Mars
exploration elaborate on its seasonal dust storms that may cover the whole red
planet. In fact, such massive dust storms regularly occur on Mars (Andersson
et al, 2015) and are becoming a major concern of academic research in the real
world. For example, Battalio and Wang (2021) conducted a comprehensive study
on the distribution and evolution of dust storms over 8 Mars years. They found
that non-global dust events frequently occur during
$L_{s}=140^{\circ}-250^{\circ}$ in the primary season as well as
$L_{s}=300^{\circ}-360^{\circ}$ (northern hemisphere) and
$L_{s}=10^{\circ}-70^{\circ}$ (southern hemisphere) in the second seasons.

Figure 1: Dusty and clean images of Mars captured by the MoRIC of Tianwen-1
orbiter. a Dusty images with heavy dust storms. b. Clean images with weak or
no dust storms.

The same phenomena have been observed by the latest Martian orbiter Tianwen-1
(see Fig. 1), which started to capture remote sensing images of Mars with its
Moderate Resolution Imaging Camera (MoRIC) in Nov. 2021 and plans to collect
color images of all areas of Mars with the resolution of $\sim$100 meters
before the end of Jun. 2022. These images will be used to construct a digital
global topographic map of Mars, while the construction process is now severely
delayed by the frequently occurring dust storms. Although obvious dust storms
can be found only in less than 5% of existing images captured by MoRIC, they
remarkably decrease the image contrast and blur/occlude features around edges,
leading to “holes” in the global topographic map. If such dust storms can be
visually eliminated or alleviated from the remote sensing images, the progress
of mapping the topography and geomorphology of Mars can be greatly boosted (Li
et al, 2021a), and the completeness of the global topographic map of Mars can
be enhanced as well. Many other subsequent tasks, such as landing site
selection (Liu et al, 2022), mineral exploration (Ehlmann et al, 2008a, b) and
water/ice activity discovery (Bibring et al, 2004; Bandfield, 2007), can also
benefit from cleaner remote sensing images.

Unfortunately, the removal of Martian dust storms from remote sensing images
is still a challenging task, since many core mechanisms and key parameters
about the formation, variation, and composition of Martian dust storms are
still far from fully explored. Although existing studies of Martian dust
storms have proposed many observations and assumptions about the atmospheric
dynamics (Vandaele et al, 2019), tidal transport of dust (Wu et al, 2020),
climate evolution (Leovy, 2001), water loss to space (Chaffin et al, 2021),
and hydrogen escape (Heavens et al, 2018), none of them attacks the problem of
Martian dust storms from a “visual” perspective:

  1. 1.

How to formulate the influence of dust storms in the imaging process?

  2. 2.

Which key parameters regularize the imaging quality degradation?

  3. 3.

How to remove the dust storms and roll back the degradation process?

## Problem Formulation

To address these challenges, we refer to a similar imaging degradation
phenomenon on the Earth (see Fig. 2a). In hazy weathers, Rayleigh scattering
occurs around small particles floating in the atmosphere (e.g., hydrone and
aerogel molecules) (McCartney, 1976). After the light transmission, the
imaging degradation process can be formulated as in (Li et al, 2021b):

$H(x,\lambda)=C(x,\lambda)\cdot{}T(x,\lambda)+L(\lambda)\cdot{}(1-T(x,\lambda)),$
(1)

where $C(x,\lambda)$ and $H(x,\lambda)$ are the clean and hazy visual signals
at location $x$ and wavelength $\lambda$, respectively. $L(\lambda)$ indicates
the global atmospheric light at wavelength $\lambda$. At a pixel $x$,
$T(x,\lambda)=\exp({-\beta(\lambda)\cdot{}S(x)})$ denotes the light
transmission ratio, in which $\beta(\lambda)$ is a positive scattering weight
at wavelength $\lambda$ and $S(x)$ is the scene depth. From Eq. 1, we can see
that a hazy image is formed as a depth-related blending of the clean image and
the global atmospheric light. At zero depth we have $T(x,\lambda)=1$ and
$H(x,\lambda)=C(x,\lambda)$. At infinite depth we have $T(x,\lambda)=0$, and
$H(x,\lambda)$ is solely determined by the global atmospheric light
$L(\lambda)$.

Figure 2: Image degradation process. a. Hazy weather on Earth; b. Dusty
weather on Mars.

With Eq. 1 in mind, we turn to the dusty weather on the Mars (see Fig. 2b).
Considering that dust storms may be non-uniformly distributed in the whole
view of the orbiter (see Fig. 1a), and the multi-grained iron oxides suspended
in the Martian atmosphere mainly have Mie scattering processes (Gao et al,
2022), we rewrite Eq. 1 to describe the Martian image degradation process as

$H(x,\lambda)=C(x,\lambda)\cdot{}\hat{T}(x,\lambda)+\hat{L}(\lambda)\cdot{}(1-\hat{T}(x,\lambda)),$
(2)

where $\hat{T}(x,\lambda)=\exp(-\hat{\beta}(\lambda)\cdot{}D(x))$ denotes the
light transmission ratio at the location $x$ when the light with wavelength
$\lambda$ travels through a local dust storm with the density $D(x)$ and the
positive scattering weight $\hat{\beta}(\lambda)$. $\hat{L}(\lambda)$ is also
the global atmospheric light at wavelength $\lambda$. Its multiplication with
$1-\hat{T}(x,\lambda)$ denotes the scattered atmosphere light and can be
intuitively viewed as the sun light reflected by the dust storm. From Eq. 2,
we can see that $\hat{T}(x,\lambda)=1$ and $H(x,\lambda)=C(x,\lambda)$ when no
dust storm exists, and $\hat{T}(x,\lambda)\approx{}0$ and
$H(x,\lambda)\approx{\hat{L}(\lambda)}$ when heavy dust storm occurs.

By formulating the influence of dust storms in the imaging process with Eq. 2,
we can see that the whole process is controlled by two key parameters,
including $\hat{T}(x,\lambda)$ and $\hat{L}(\lambda)$. Unfortunately, these
key parameters are still unknown, which prevent the direct inference from the
degraded image $H(x,\lambda)$ to the desired clean image $C(x,\lambda)$. To
address this issue, we propose a deep learning method to resolve the problem
in a data-driven manner.

## Method

To remove Martian dust storms with deep learning, the first key step is to
synthesize dusty images from clean ones to supervise the learning process. In
particular, the data synthesis process must be fully controllable without
changing any other image content except adding dust storms, which prevent us
from using unsupervised techniques such as Cycle-Consistent Adversarial
Networks (Zhu et al, 2017). Therefore, we refer to Eq. 2 and synthesize dusty
images with random parameters.

Figure 3: Comparison of distributions of different types of Mars images. a.
realistic clean and dusty images have remarkably different feature
distributions. b. Synthetic and realistic dusty images have similar feature
distributions.

First, we collect 6,746 Martian images captured by the MoRIC of Tianwen-1
between Nov. 2021 and Mar. 2022. Each image has a resolution of
$3072\times{}4096$. From these images, we manually select 500 images with dust
storms (about 7.4%) and 500 clean images (see Fig. 1b). After that, we train a
deep classification model with the backbone ResNet-50 (He et al, 2016) to
separate clean and dusty images. Based on the extracted features, we visualize
the data distribution in Fig. 3a by using the t-SNE algorithm (Van der Maaten
and Hinton, 2008), in which clean and dusty images have two different
distributions. Thus the objective of data synthesis can be described as
generating synthetic dusty images using Eq. 2 that have the same distribution
as realistic dusty images.

Toward this end, we propose to directly generate random transmission ratio
maps $\hat{T}(x,\lambda)\in{}[0,1]$ that inherently cover the influences of
both $\hat{\beta}(\lambda)$ and $D(x)$ by using the 2D Perlin noise (Perlin,
1985). As shown in Fig. 4, the Perlin noise have similar visual effects as
dust storms on the Mars. It is controlled by four parameters, including Scale
(the bounds of noise), Octaves (level of noise details), Lacunarity (levels of
detail bias in each octave), and Persistence (influence of each octave to the
overall shape). For each clean image, we randomly generate seven Perlin noise
maps and re-weight each map with random
$\alpha\in\\{0.4,0.5,0.6,0.7,0.8,0.9,1.0\\}$. Let $M(x)\in{}[0,1]$ be the
Perlin noise at $x$, we generate $\hat{T}(x,\lambda)\in{}[0,1]$ as

$\hat{T}(x,\lambda)=1-\alpha\cdot{}M(x).$ (3)

In addition, the global atmospheric light $\hat{L}(\lambda)$ can be
intuitively viewed as the sunlight reflected by the dust storm with infinite
density. As in (Tan, 2008), the intensity of sun light at any wavelength can
be roughly estimated as the maximum of the clean image $C(x,\lambda)$, and
$\hat{L}(\lambda)$ further incorporated the reflexivity of the dust storm,
denoted as $\phi(\lambda)$, into account:

$L(\lambda)=\phi(\lambda)\cdot{}\max_{x,\lambda_{0}}{}C(x,\lambda_{0}).$ (4)

In practice, we manually select image patches covered by heavy dust storms
from the 500 realistic dusty images, denoted as $\mathbb{H}_{0}$. As a result,
the reflexivity can be heuristically estimated as

$\phi(\lambda)=\frac{1}{\|\mathbb{H}_{0}\|}\sum\limits_{H\in{}\mathbb{H}_{0}}\left(\frac{1}{\|H\|}\sum\limits_{x\in{}H}\frac{H(x,\lambda)}{\max\limits_{\lambda{}_{0}}H(x,\lambda_{0})}\right),$
(5)

where $\|H\|$ denotes the number of pixels in the image patch $H$, and
$\|\mathbb{H}_{0}\|$ denotes the number of selected patches covered by heavy
dust storms.

Figure 4: Examples of Perlin noise maps generated with different parameters.

By incorporating Eq. 3, Eq. 4 and Eq. 5 into Eq. 2, we can generate 3500
synthetic dusty images from the 500 clean images with the roughly estimated
values of key parameters and randomly generated Perlin noise maps. For these
synthetic images, we extract their deep features by using the same ResNet-50
model that can distinctively separate clean and dusty images. As shown in Fig.
3b, these synthetic dusty images have similar distributions with realistic
ones and thus forming 3500 clean-dusty image pairs.

Given these image pairs, we can train a deep network with the encoder-decoder
architecture (see Fig.5). In this network, image features are extracted by
three dense depth-wise separable convolution (DDSC) modules after the image is
down-sampled twice by a factor of 2 to gradually reduce spatial resolution and
extend feature channels. Note that each DDSC module contains multiple depth-
wise separable convolutions cascaded in a densely connected manner for
efficient feature extraction. Besides, we adopt the FA block from (Qin et al,
2020) to re-weight different channels and different spatial locations of
feature maps. Finally, two up-sampling layers and one convolution layer are
used to generate a dust-free image that has the same size with the input
image. In this manner, the reverse image degradation process of Eq. 2 is
inherently modeled by the network architecture.

Figure 5: The architecture of the proposed deep neural networks.

To train this network, we randomly crop 512 $\times$ 512 patches from the
original image as the input. Data augmentation is also used by randomly
rotating and horizontally flipping the patches. Other training details
include: using the AdamW optimizer with a batch size of 8 and a learning rate
of 0.0001 and training the network with 180 epochs.

## Results

To validate the effectiveness of the trained deep model, we test it on the 500
images captured by Tianwen-1 with heavy dust storms. The generated dust-
removal images can be found in Fig. 6, from which we find that the trained
network can effectively remove dust storms without changing any other content
of the images. Even when severe dust storms occur, the model still outputs
clear enhancement results. This allows the orbiter to see the Martian
landscape covered by dust storms more clearly and reveal more topographical
and geomorphological details of Mars.

Figure 6: Representative results generated by the proposed approach. a,c.
Images captured by Tianwen-1 with heavy dust storms. b,d. Dust-removal
results.

Beyond qualitative results, we also choose a no-reference image quality
assessment metrics proposed by (Choi et al, 2015) to quantitatively evaluate
the image quality before and after the dust-removal operation. The metric,
denoted as FADE, measures the intensity of haze in the image, and we use it to
assess the intensity of dust storms. The smaller a FADE score, the less
haze/dust. We find that the average FADE score of the 500 realistic clean
images is only 0.5022, and the average FADE scores of the 500 realistic dusty
images and the 3500 synthetic dusty images are 0.7633 and 0.7627,
respectively. These results indicate that the metric FADE is suitable to
assess the intensity of dust storms, and the synthetic images generated using
our method have almost the same distribution of dust storms as in the real-
world. Finally, the 500 dust-removal images generated by the proposed approach
have an average FADE score of 0.5907, implying that the trained network has an
impressive performance in removing dust storms from Martian images. These
results also shows the powerful capability of the deep models in extracting
dust-irrelevant features. By using three consecutive DDSC modules, robust
features are extracted to encode the content of the Martian landscape, while
dust storms are gradually filtered out as some kinds of “noise.” As a result,
images recovered from these features contain only the topographic details of
the landscape without reconstructing such “noise,” leading to a dust-free
decoding results.

## Discussion

Based on the results, the key issue needs to be further discussed is why the
deep learning approach performs impressively in resolving the dust-removal
task. By looking into the convolutional mechanisms of deep models and the
image degradation process of dust storms, we find two possible explanations.

First, the influence of dust storms at each pixel acts as a random noise.
Since the density of dust storms rarely reaches the theoretical maximum (i.e.,
infinity), every pixel is actually a mixture of ground-reflected and
atmosphere-scattered light with a random non-zero weight. However, the ground-
reflected light intensities of nearby pixels are tightly correlated since
Martian landscapes are highly structured, which enables the convolution
kernels to remove the random noise from atmosphere-scattered light after
addition and multiplication operations and to maintain the local structure.

Second, the training process of the deep model is fully supervised at each
pixel. Inspired by the haze formation process on Earth, we formally formulate
the image degradation process when dust storms occur. In this manner, massive
synthetic data can be generated with near-realistic image contents. In this
manner, the formation and removal of dust storms are actually regularized by
this mathematical equation, while such mathematical regularization is
inherently remembered by the deep models after the training process. That also
explains why the trained deep model can successfully remove the dust storms
without changing any other image content.

To sum up, the success of the proposed approach mainly comes from the nature
of convolutional operations and the mathematical problem formulation, which
resolves the three challenges in eliminating Martian dust storms from a visual
perspective and roll back the degradation process. In the future work, we will
further incorporate the human-in-the-loop mechanism to enable the dust-removal
model eliminate various types of dust storms more accurately by actively
learning from human interactions and annotations. Such a dust-removal
technique will also be released as a tool to boost the exploration of Mars
surface structure and mapping its global morphology to form a digital global
topographic map of Mars using Tianwen-1.

Acknowledgments

This work was supported by the National Natural Science Foundation of China
(No. 62132002, 61922006). The Tianwen-1 data used in this work is processed
and produced by “Ground Research and Application System (GRAS) of China’s
Lunar and Planetary Exploration Program, provided by China National Space
Administration (http://moon.bao.ac.cn)”.

## References

  *   * Andersson et al (2015) Andersson L, Weber T, Malaspina D, et al (2015) Dust observations at orbital altitudes surrounding mars. Science 350(6261):aad0398
  * Bandfield (2007) Bandfield JL (2007) High-resolution subsurface water-ice distributions on mars. Nature 447(7140):64–67
  * Battalio and Wang (2021) Battalio M, Wang H (2021) The mars dust activity database (mdad): A comprehensive statistical study of dust storm sequences. Icarus 354:114,059
  * Bibring et al (2004) Bibring JP, Langevin Y, Poulet F, et al (2004) Perennial water ice identified in the south polar cap of mars. Nature 428(6983):627–630
  * Chaffin et al (2021) Chaffin M, Kass D, Aoki S, et al (2021) Martian water loss to space enhanced by regional dust storms. Nature Astronomy 5(10):1036–1042
  * Choi et al (2015) Choi LK, You J, Bovik AC (2015) Referenceless prediction of perceptual fog density and perceptual image defogging. IEEE Transactions on Image Processing 24(11):3888–3901
  * Ehlmann et al (2008a) Ehlmann BL, Mustard JF, Fassett CI, et al (2008a) Clay minerals in delta deposits and organic preservation potential on mars. Nature Geoscience 1(6):355–358
  * Ehlmann et al (2008b) Ehlmann BL, Mustard JF, Murchie SL, et al (2008b) Orbital identification of carbonate-bearing rocks on mars. Science 322(5909):1828–1832
  * Gao et al (2022) Gao X, Xie L, Dou X, et al (2022) Effects of charged martian dust on martian atmosphere remote sensing. Remote Sensing 14(9):2072
  * He et al (2016) He K, Zhang X, Ren S, et al (2016) Deep residual learning for image recognition. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp 770–778
  * Heavens et al (2018) Heavens NG, Kleinböhl A, Chaffin MS, et al (2018) Hydrogen escape from mars enhanced by deep convection in dust storms. Nature Astronomy 2(2):126–132
  * Leovy (2001) Leovy C (2001) Weather and climate on mars. Nature 412(6843):245–249
  * Li et al (2021a) Li C, Zhang R, Yu D, et al (2021a) China’s mars exploration mission and science investigation. Space Science Reviews 217(4):1–24
  * Li et al (2021b) Li H, Li J, Zhao D, et al (2021b) Dehazeflow: Multi-scale conditional flow network for single image dehazing. In: Proceedings of the 29th ACM International Conference on Multimedia, pp 2577–2585
  * Liu et al (2022) Liu J, Li C, Zhang R, et al (2022) Geomorphic contexts and science focus of the zhurong landing site on mars. Nature Astronomy 6(1):65–71
  * Van der Maaten and Hinton (2008) Van der Maaten L, Hinton G (2008) Visualizing data using t-sne. Journal of machine learning research 9(11)
  * McCartney (1976) McCartney EJ (1976) Optics of the atmosphere: scattering by molecules and particles. New York
  * Perlin (1985) Perlin K (1985) An image synthesizer. ACM Siggraph Computer Graphics 19(3):287–296
  * Qin et al (2020) Qin X, Wang Z, Bai Y, et al (2020) Ffa-net: Feature fusion attention network for single image dehazing. In: Proceedings of the AAAI Conference on Artificial Intelligence, pp 11,908–11,915
  * Robinson (1993) Robinson KS (1993) Red Mars. Spectra, Waukegan
  * Scott (2015) Scott R (2015) The martian. 20th Century Studios
  * Tan (2008) Tan RT (2008) Visibility in bad weather from a single image. In: 2008 IEEE Conference on Computer Vision and Pattern Recognition, IEEE, pp 1–8
  * Vandaele et al (2019) Vandaele AC, Korablev O, Daerden F, et al (2019) Martian dust storm impact on atmospheric h2o and d/h observed by exomars trace gas orbiter. Nature 568(7753):521–525
  * Wu et al (2020) Wu Z, Li T, Zhang X, et al (2020) Dust tides and rapid meridional motions in the martian atmosphere during major dust storms. Nature Communications 11(1):1–10
  * Zhu et al (2017) Zhu JY, Park T, Isola P, et al (2017) Unpaired image-to-image translation using cycle-consistent adversarial networks. In: Proceedings of the IEEE International Conference on Computer Vision, pp 2223–2232