# Sym-Q: Adaptive Symbolic Regression via Sequential Decision-Making

Yuan Tian  Wenqi Zhou  Hao Dong  David S. Kammer  Olga Fink

###### Abstract

Symbolic regression holds great potential for uncovering underlying
mathematical and physical relationships from empirical data. While existing
transformer-based models have recently achieved significant success in this
domain, they face challenges in terms of generalizability and adaptability.
Typically, in cases where the output expressions do not adequately fit
experimental data, the models lack efficient mechanisms to adapt or modify the
expression. This inflexibility hinders their application in real-world
scenarios, particularly in discovering unknown physical or biological
relationships. Inspired by how human experts refine and adapt expressions, we
introduce Symbolic Q-network (Sym-Q), a novel reinforcement learning-based
model that redefines symbolic regression as a sequential decision-making task.
Sym-Q leverages supervised demonstrations and refines expressions based on
reward signals indicating the quality of fitting precision. Its distinctive
ability to manage the complexity of expression trees and perform precise step-
wise updates significantly enhances flexibility and efficiency. Our results
demonstrate that Sym-Q excels not only in recovering underlying mathematical
structures but also uniquely learns to efficiently refine the output
expression based on reward signals, thereby discovering underlying
expressions. Sym-Q paves the way for more intuitive and impactful discoveries
in physical science, marking a substantial advancement in the field of
symbolic regression.

Symbolic Regression, Reinforcement Learning

##  1 Introduction

Symbolic regression is a powerful form of regression analysis that searches
the space of mathematical expressions to find the expression that best fits a
given dataset. Unlike traditional regression models that fit data to pre-
specified equations, symbolic regression is capable of discovering the
underlying equations or relationships between variables. This capability can
lead to a deeper understanding of the inherent structure and dynamics of the
data and is particularly important in fields where the relationships between
variables are complex and not well understood, as it provides a tool to
uncover the form of the relationship without prior assumptions. Recently,
symbolic regression has been instrumental in uncovering new relationships,
such as astrophysical scaling relations (Wadekar et al., 2022) and analytical
models of exoplanet transit spectroscopy (Matchev et al., 2022). However, a
significant challenge in symbolic regression is its inherent complexity. This
complexity grows rapidly with the length of the symbolic expressions, making
it a computationally demanding NP-hard problem (Li et al., 2022).

Figure 1: Overview of the proposed framework. A. The expression and its
corresponding expression tree. B. The proposed Sym-Q agent supports both
offline training with ground truth human knowledge and online searching with
reward signals. $\tau_{circle}$ represents the agent trajectories within the
symbolic environment. C. The Sym-Q architecture and step-wise decision-making
mechanism.

Symbolic regression has primarily been explored through two approaches: online
search methods like Genetic programming (GP) and transformer-based models
trained with extensive data. GP, a leading online search method (Forrest,
1993; Schmidt & Lipson, 2009; Augusto & Barbosa, 2000; Gustafson et al., 2005;
Bladek & Krawiec, 2019), iteratively evolves each generation of mathematical
expressions to match the observation points through selection, crossover, and
mutation. Besides, previous reinforcement learning-based (RL) methods, like
deep symbolic optimization (DSO) (Petersen et al., 2019), optimize expressions
based on reward signals during online search. Some recent methods even combine
GP with RL, leveraging the advantages of both approaches (Mundhenk et al.,
2021). However, despite their effectiveness, these online search methods treat
symbolic regression as an instance-based problem, requiring model training or
searching from scratch for each case, leading to inefficiency and high
computational demands.

Building upon the advancements facilitated by large-scale language models
(Vaswani et al., 2017), recent efforts in symbolic regression have focused on
using transformer-based models, which are trained on extensive datasets (Li et
al., 2022; Vastl et al., 2022; Kamienny et al., 2022; Bendinelli et al.,
2023). Although transformer-based models have recently achieved significant
success in symbolic regression, showing considerable capability in finding the
underlying equations for many cases, they still face challenges in
generalizability and adaptability. These models lack efficiency in adapting
their output expressions after training. This lack of flexibility severely
constrains their utility in practical applications, as it requires the
training dataset to comprehensively cover relevant expressions. Such rigidity
is particularly limiting in fields like physics or biology, where the ability
to uncover new relationships hinges on the model’s flexibility and
adaptability.

Transformer-based models in symbolic regression can sometimes produce near-
correct expressions with minor errors, like a single incorrect operator.
Rather than re-optimizing the entire sequence from beginning, an efficient
model be able to pinpoint and correct just the erroneous element.

To address this challenge, we introduce a novel approach by conceptualizing
symbolic regression as a Markov decision process (MDP) for the first time. In
our new framework, the model assesses the current expression tree and observed
data at each step to choose an action that maximizes cumulative rewards like
the coefficient of determination $(R^{2})$ (Glantz et al., 2001). This
reformulation transforms the large action space ($m^{l}$, where $m$ is the
number of total possible operations for each step and $l$ is the maximum
decision depth) into a more manageable scope, reducing the actions per step to
just the number of total possible operations $m$. Consequently, this novel
formulation facilitates efficient, step-wise updating and adapting of the
expression (Tian et al., 2020).

Under this formulation, in this work, we introduce Symbolic Q-network (Sym-Q),
a novel RL-based method for symbolic regression as illustrated in Figure 1.
This approach, under the MDP formulation, innovatively combines online RL
techniques with pretrained transformer-based models. Sym-Q is the first
approach to employ sequential decision-making for generating and refining
expression trees, representing a pioneering step in the intersection of RL and
symbolic regression. It leverages the strengths of both RL and advanced
transformer models to iteratively explore and optimize the space of
mathematical expressions. The key contributions of our study are threefold:

  * •

We reformulate symbolic regression as a sequential decision-making process,
demonstrating its effectiveness.

  * •

Building upon this new problem formulation, we propose the Sym-Q model. Sym-Q
is trained via offline RL and is also compatible with online RL methods,
enabling efficient refinement and adaptation of output expressions in response
to reward signals.

  * •

Empirically, Sym-Q outperforms current state-of-the-art transformer-based
models, particularly in recovering expression structures and refining output
expressions based on rewards. We provide an in-depth analysis of the model’s
performance, offering insights for future training data design and method
improvements. Our source code will be publicly available.

##  2 Related works

The field of symbolic regression has been primarily explored in two
methodological trajectories: online searching methods, and transformer-based
models. Genetic programming-based approaches are the most widely used in the
online searching category and include many notable works such as (Forrest,
1993; Schmidt & Lipson, 2009; Augusto & Barbosa, 2000; Gustafson et al., 2005;
Bladek & Krawiec, 2019). By far the most popular commercial software Eureqa
(Dubčáková, 2011) is the most successful application based on GP methods.

Within the field of online search methods for symbolic regression, RL-based
approaches represent a significant advancement (Petersen et al., 2019; Tenachi
et al., 2023), as highlighted by Petersen et al. (2019) in their introduction
of deep symbolic optimization (DSO). DSO employs a risk-seeking policy
gradient to effectively train a recurrent neural network (RNN), which then
generates a probability distribution over the space of mathematical
expressions. Additionally, this field has witnessed the advent of diverse
search methodologies, such as normalizing flows (Udrescu & Tegmark, 2020) and
the combination of genetic programming with RL (Mundhenk et al., 2021), aiming
to capitalize on the strengths of both. Despite their potential, these methods
typically require searches or training for each instance from scratch, leading
to computational challenges and stochastic variability.

The second methodological approach in symbolic regression focuses on
transformer-based models trained on extensive datasets and has shown promise
in identifying underlying equations in various scenarios. Recent works,
including those by Li et al. (2022); Vastl et al. (2022); Kamienny et al.
(2022); Bendinelli et al. (2023), demonstrate the effectiveness of transformer
models in this field. Models like NeSymReS (Biggio et al., 2021) and
SymbolicGPT have pioneered the use of Set transformers (Devlin et al., 2018;
Lee et al., 2019) and GPT frameworks (Valipour et al., 2021), respectively,
for encoding observed data points and generating expressions. Li et al. (2022)
introduce a novel approach employing PointMLP (Ma et al., 2022) for local and
global feature extraction and a joint learning mechanism that combines
supervised contrastive learning with the guidance of the preorder traversal of
expression trees. Neural Symbolic Regression with Hypothesis (NSRwH)
(Bendinelli et al., 2023) proposes the inclusion of assumptions about the
ground truth expression structure in the prediction process. While Kamienny et
al. (2022) proposes to directly output complete expressions, including
constants, these mentioned models’ generally face limitations in adapting
output expressions efficiently, which restricts their practical application in
real-world scenarios. Notably, existing learning-based works, whether RL-based
or transformer-based, do not facilitate step-wise updates and refinements, a
gap that our approach aims to address.

##  3 Method

###  3.1 Overview

In this work, we propose Sym-Q, an RL-based method for addressing large-scale
symbolic regression problems. The overview of the proposed Sym-Q is depicted
in Figure 1. By reformulating symbolic regression as an MDP, the agent is
trained to understand the observations and sequentially expand and refine the
expression tree. As illustrated in Figure 1.B, Sym-Q is capable of learning
from both optimal demonstrations, and sub-optimal reward signals. Initially,
the agent is trained using the proposed offline RL method and can further
refine its policy through our proposed off-policy RL method during online
searching.

Figure 1.C details the decision-making process. At each time step, the agent
observes the encoded data points and the encoded preceding expression tree.
Utilizing this encoded information, the agent identifies the placeholder node
(shown as red node in the tree) that needs to be replaced. Then the agent
sequentially replaces the placeholders and expands the expression tree until
termination. This process is performed $128$ times using the beam search
strategy, resulting in $128$ distinct skeletons.

The final expression is then derived by applying the
Broyden–Fletcher–Goldfarb–Shanno (BFGS) (Fletcher, 2000) optimization,
focusing on the constant terms in the expressions. For detailed
implementation, the pseudo-code of our deployment strategy is available in
Appendix A.

###  3.2 Expression tree

An expression tree, as described by Hopcroft et al. (2006), is a structure
where internal nodes symbolize mathematical operations, such as addition,
subtraction, and logarithm, while leaf nodes represent constants or variables.
Due to its tree structure, it is particularly suitable for performing
sequential decision-making. This structure, visually depicted in Figure 1.A,
facilitates the systematic assembly and analysis of mathematical expressions,
especially when their exact forms are initially unknown.

###  3.3 Symbolic regression as a sequential decision-making task

Symbolic regression involves searching for a sequence
$a_{1:T}=a_{1},a_{2},...,a_{T}$ and constructing it into an expression
represented by an expression tree to fit observational data. In this context,
$a_{t}$ represents an operation within a specified search space at time step
$t$.

Previous large-scale transformer-based models approach the symbolic regression
problem as a machine translation problem. They use feature extractors on the
input data points and then decode the features into expression skeletons in an
auto-regressive fashion, which means they produce the entire symbolic tree in
a single output. In our study, we reinterpret the task of symbolic regression
as an Markov decision process (MDP). An MDP is defined by the tuple
$(S,A,r,P,\rho)$, where $S$ denotes the set of states that accurately
represent the current situation, $A$ is the set of possible actions, $r(s,a)$
is the reward function, $P(s^{\prime}|s,a)$ describes the transition
probabilities, and $\rho(s)$ is the initial state distribution.

At each time step $t$, the model selects an action $a_{t}$ based on the policy
$\pi_{\theta}(a_{t}|s_{t})$, with $s_{t}$ being the current state. This state
encapsulates both the observed data points $({X},{y})$ and the encoded
structure of the current tree, derived from previous actions $a_{1:t-1}$ . The
action $a_{t}$ is chosen to determine the next operation in the sequence. The
symbolic expression $a_{1:T}=a_{1},a_{2},...,a_{T}$ is then constructed as a
trajectory of the MDP, represented as
$\tau={s_{1},a_{1},s_{2},a_{2},...,s_{T},a_{T}}$. In our model, rewards are
granted only at the termination state, with all other states receiving a
reward of $0$. This approach is consistent with strategies used in other
domains such as various games (Silver et al., 2017, 2016), and robotics
(Vecerik et al., 2017; Riedmiller et al., 2018). As a reward measure, we
utilize the coefficient of determination $(R^{2})$ (Glantz et al., 2001),
although it is not the only option for this task. The $R^{2}$ is defined as:

$R^{2}=1-\frac{\sum_{i=1}^{n}(y_{i}-\hat{y}_{i})^{2}}{\sum_{i=1}^{n}(y_{i}-\bar{y})^{2}},$
(1)

where $y_{i}$ and $\hat{y}_{i}$ represent the ground-truth and predicted
values for point $i$, respectively. $\bar{y}$ is the average of the $y_{i}$
values across all data points, and $n$ denotes the number of observed data
points.

Under the introduced MDP formulation, the goal of an RL-based agent is to find
a policy $\pi$ that maximizes the cumulative reward $J(\pi)$:

$J(\pi)=\mathbb{E}_{\tau\sim\rho_{\pi}}{\sum_{t=0}^{\infty}r(s_{t},a_{t})}.$
(2)

Online RL algorithms require extensive trial and error to decipher the rules
of symbolic regression. This task is particularly challenging as it targets
only optimal solutions, with all sub-optimal solutions receiving significantly
lower rewards. The sparsity of the optimal reward signal exacerbates the
inefficiency of searching. Offline RL offers a learning paradigm that merges
the ability of supervised learning to leverage existing data with RL’s
capacity to optimize arbitrary rewards and leverage temporal compositionality
(Levine et al., 2020; Kostrikov et al., 2021; Kumar et al., 2020; Janner et
al., 2021). For offline RL, while the ultimate objective remains the
optimization of the function specified in Equation 2, there is a fundamental
change: the agent is unable to interact directly with the environment or
collect new transitions through the behavior policy. Instead, the learning
process relies on a predetermined, static dataset of transitions denoted as
$D={(s^{i}_{t},a^{i}_{t},s^{i}_{t+1},r^{i}_{t})}$. This dataset replaces real-
time interaction, serving as the primary source for the learning algorithm.

###  3.4 Symbolic Q-network

Sym-Q is specifically designed to enable sequential decision-making and
efficient step-wise updates with online and offline data. It comprises three
integral modules:

  1. 1.

Point set encoder $E^{points}_{\phi}({{X,y}})$ which transforms the point sets
${(X,y)}=((x_{1}^{1},x_{1}^{2},y_{1}),(x^{1}_{2},x_{2}^{2},y_{2}),...,(x^{1}_{n},x_{n}^{2},y_{n}))$
associated with each equation into a latent space, resulting in a fixed-size
latent representation denoted as $z^{p}\in\mathbb{R}^{1\times K_{p}}$, where
$K_{p}$ indicates the dimension of the latent variable;

  2. 2.

Expression tree encoder $E^{tree}_{\psi}(M_{t})$ which maps the current tree’s
one-hot encoded matrix $M_{t}$ into another fixed-size latent representation
$z^{t}\in\mathbb{R}^{1\times K_{t}}$, where $K_{t}$ indicates the dimension of
the variable;

  3. 3.

Q-network $Q_{\theta}(s_{t})$ which is responsible for calculating the Q-value
for each potential operation, given the combined latent state
$s_{t}=(z^{p},z^{t})$.

This structure enables Sym-Q to approach symbolic regression in a novel and
efficient way. The overall framework of Sym-Q is depicted in Figure 1. For
further details of the neural network architectures and hyperparameters used,
please refer to Appendix B.

###  3.5 Efficient conservative offline Q-learning for symbolic regression

In training the Sym-Q through offline RL, a significant challenges is the
overestimation of values due to distributional shifts between the dataset and
the learned policy. To address this, we adopt Conservative Q-learning (CQL)
(Kumar et al., 2020), which minimizes the values of state-action pairs outside
the training data distribution, while simultaneously maximizing the values
within this distribution. Given the sparse reward structure of symbolic
regression, where rewards are typically given only at the completion of each
trajectory and most sub-optimal trajectories receive lower rewards, we have
adapted CQL to this context. Our modified version of CQL is customized for
symbolic regression, aiming to effectively leverage optimal offline data in
this domain:

$J(\theta,\psi)=-\mathbb{E}_{(s,a^{i})\sim\mathcal{D}}\left[\log\left(\frac{e^{Q_{\theta,\psi}(s,a^{i})}}{\sum_{j=1}^{m}e^{Q_{\theta,\psi}(s,a^{j})}}\right)\right],$
(3)

where $Q_{\theta,\psi}$ represents the network and its weight parameters,
including those beyond the points encoder. $D$ represents the training
dataset, $a^{i}$ is the action demonstrated at state $s$, and ${a^{j}|j\neq
i,1\leq j\leq m}$ are the non-demonstrated actions, regarded as suboptimal.

Building on the relationship between conservative Q-learning loss and cross-
entropy loss, our proposed objective function represents the log probability
of the optimal action. This mirrors the emphasis of cross-entropy on
identifying the correct class. The inclusion of the softmax function in the
denominator serves to normalize the Q values across all possible actions,
effectively converting them into a probability distribution.

By adopting this objective function, our approach more closely aligns with the
principles of traditional supervised learning. This alignment could make the
approach more intuitive and straightforward to implement, particularly in
scenarios where the correct action at each state is well-defined, such as in
symbolic regression. A detailed explanation of the algorithm can be found in
Appendix C.

###  3.6 Supervised contrastive learning for points encoder

Following a similar idea to Li et al. (2022), we implemented the supervised
contrastive loss for expressions and data points that share the same skeleton.
This approach uses the skeleton of expressions as category labels to enrich
supervisory information. More specifically, expressions and their
corresponding latent points encoding $z^{p}_{i}$, which belong to the same
skeleton, are grouped together in the embedding space. At the same time, we
ensure that the latent points encoding $z^{p}_{j}$ from different skeletons
are distanced from each other. The supervised contrastive loss is defined as:

$\displaystyle J(\phi)$
$\displaystyle=\sum_{i=1}^{N}\frac{-1}{|\mathcal{P}(i)|}\sum_{p\in\mathcal{P}(i)}\log$
(4) $\displaystyle\frac{\exp(E^{points}_{\phi}({X,y})^{i})\cdot
E^{points}_{\phi}({X,y})^{p}/\tau)}{\sum_{j=1}^{N}1_{[j\neq
i]}\exp(E^{points}_{\phi}({X,y})^{i}\cdot
E^{points}_{\phi}({X,y})^{j}/\tau)},$

where $\mathcal{P}(i)$ represents the set of indices for all positives in the
multiviewed batch distinct from $i$, and $|\mathcal{P}(i)|$ is its
cardinality; $\tau$ is an adjustable scalar temperature parameter controlling
class separation; $N$ represents the mini-batch size;
$z^{p}_{i}=E^{points}_{\phi}({X,y})^{i}$ is the points embedding of the sample
$i$ in a batch, $z^{p}_{p}=E^{points}_{\phi}({X,y})^{p}$ is the points
embedding of a positive sample sharing the same skeleton as sample $i$. The
overall loss is then given by:

$\displaystyle J(\phi,\theta,\psi)=J(\theta,\psi)+\alpha J(\phi),$ (5)

where $\alpha$ is a scalar weighting hyperparameter tuned during the training
stage. More details about the parameters settings and network architectures
can be found in Appendix B.

###  3.7 Online searching and updating

Acknowledging the limitations in current approaches, it becomes evident that a
well-generalized model capable of discovering any possible equation would
require an impractically large training dataset, making it infeasible for
real-world applications. Therefore, a more practical solution lies in
developing a pretrained model that is able to conduct efficient and effective
online searches when faced with new formulas. In line with this idea, we
propose to build on the offline reinforcement learning objective function
detailed in Equation 5 and introduce an off-policy risk-seeking REINFORCE
method. This method, inspired by DSO (Petersen et al., 2019), focuses on
maximizing best-case performance and setting aside less optimal scenarios.
Such a strategy is not only pragmatic but has also been proven effective in
enhancing the adaptability and responsiveness of models in the realm of
symbolic regression:

$\displaystyle J(\theta)=$ $\displaystyle-\mathbb{E}_{(s,a^{i},R_{\tau}){\sim
M}}(R_{\tau}-\beta
R^{*})\log\left(\frac{e^{Q_{\theta}(s,a^{i})}}{\sum_{j=1}^{m}e^{Q_{\theta}(s,a^{j})}}\right),$
(6)

where $R_{\tau}$ is the online cumulative reward of the corresponding Markov
tuple, $R^{*}$ is the highest cumulative reward during searching, $\beta$
controls the performance threshold of considered cases, and $M$ is the memory
buffer during online searching. Equation 5 is a special case that the
$R_{\tau}=1$. For the online searching strategy, we simply sample the action
based on normalized value distribution.

###  3.8 Optimization of the coefficients

The BFGS optimization algorithm, commonly used in previous research for
calculating coefficients within the skeleton of derived equations, is also
employed in our approach, as outlined by Li et al. (2022).

##  4 Experimental Results

###  4.1 Training datasets

For our training dataset, we created $5$ million mathematical expressions
derived from $100,000$ pre-defined expression skeletons featuring up to two
independent variables (i.e., $d=2$ and $y=f(x_{1},x_{2}))$, following the
guidelines provided by Li et al. (2022). These $100,000$ skeletons were
generated using the method described by Lample & Charton (2019), which
involves constructing a random unary-binary tree and filling its nodes with
operators and its leaves with independent variables or constants. We then
varied the constant values in each expression $50$ times and sampled $100$
random data points from each expression, adhering to the procedure outlined in
Li et al. (2022). More details about the dataset and the corresponding action
space can be found in Appendix D.

Figure 2: Training curves. These two figures demonstrate the progression of
both step accuracy (a) and equation accuracy throughout the training process
(b). The x-axis represents the training steps, while the y-axis represents
step accuracy and equation accuracy, respectively. The curves have been
smoothed using a moving average to facilitate a clearer interpretation of
trends.

###  4.2 Agent training

We trained the agent using the offline RL algorithm outlined in Section 3.
During the training, the agent demonstrated notable improvements in
performance, as indicated by the step accuracy and equation accuracy metrics,
as depicted in Figure 2. The training curves suggest continuous learning by
the agent, implying that extending the training duration could potentially
yield further enhancements in performance. We train the model on one NVIDIA
GeForce RTX 3090.

###  4.3 Evaluation on benchmark datasets

In this study, our primary focus is on generating accurate expression
skeletons, which are crucial part for deriving physical relationships. It is
important to recognize that while various combinations of operator terms can
yield a high coefficient of determination, this does not inherently assure the
precision of the final predictions or the correctness of the ultimate
expression (Li et al., 2022). We evaluate the skeleton recovery rate using the
challenging test set proposed by Li et al. (2022), which consists of the same
skeletons as in the training set but with different numerical coefficients,
named SSDNC dataset (Li et al., 2022). SSDNC dataset is currently the most
comprehensive dataset, containing 963 unseen equations. Further details about
the SSDNC dataset can be found in Appendix D. Our method is compared against
three state-of-the-art (SOTA) transformer-based supervised learning methods,
SymPointMLP (Li et al., 2022), SymbolicGPT (Valipour et al., 2021), and
NeSymReS (Biggio et al., 2021). Our method achieves a notable recovery rate of
42.7$\%$ without the beam search and an impressive 82.3$\%$ with the
implementation of the beam search strategy. As detailed in Table 1, our
proposed method surpasses the SOTA methods by up to $32.0\%$. Furthermore, it
outperforms these methods in terms of average $R^{2}$, achieving a remarkable
score of 0.95135.

We also achieve superior weighted average $R^{2}$ fitting accuracy across five
additional well-known benchmarks, including Nguyen (Uy et al., 2011),
Constant, Keijzer (Keijzer, 2003), R rationals (Krawiec & Pawlak, 2013), and
the AI-Feynman database (Udrescu & Tegmark, 2020). In addition to the three
transformer-based as outlined above, we compared with two additional online
searching approaches, DSO (Petersen et al., 2019) and the standard GP-based
symbolic regression (Koza, 1994). As detailed in the Appendix E, our method
outperforms all compared methods with a weighted average $R^{2}$ of 0.95044
across these six benchmarks. Since on these benchmarks, the skeleton recovery
rate has not been previously reported, our comparison is based solely on the
weighted average $R^{2}$. All evaluations and baseline implementation were
conducted following the setup of SymPointMLP (Li et al., 2022).

Methods | Recovery rate | $R^{2}$
---|---|---
SymbolicGPT | $50.3\%$ | 0.74087
NeSymRes | $63.4\%$ | 0.85792
SymPointML | $75.2\%$ | 0.94782
Ours | 82.3$\%$ | 0.95135
Table 1: Recovery rate of expression skeletons and $R^{2}$ on the SSDNC
benchmark. Sym-Q outperforms the state-of-the-art baselines in both skeleton
recovery rate and final fitting $R^{2}$ value. All models use the beam search
strategy with the beam size equaling 128.

##  5 Analysis of Distinct Properties of Sym-Q

Unlike previous methods that rely solely on transformer decoders, our Sym-Q
model uniquely enables step-wise analysis of its decision-making process. This
capability not only provides deeper insights into expression generation but
also offers a more accurate and detailed metric for evaluating the agent’s
performance. In the following, we will explore various aspects of our
findings: the error rate trend throughout the expression generation process,
the distribution of incorrect selections, and the specific types of errors
encountered. All analyses are conducted using the SSDNC dataset. The
computation of step-wise accuracy is based on the input of the correct
preceding expression tree encoding and the observed data point to determine if
the model’s output decision is consistent with the expected decision in the
SSDNC dataset. This is done without considering other potential solutions. The
overall step accuracy achieved is $93.5\%$.

Figure 3: The trend of error rate throughout the expression generation
process. The y-axis represents error rates, while the x-axis indicates the
progress percentage of the expression generation process, whereby $100\%$
represents a complete equation. A noticeable pattern emerges from this
representation: the middle stages of the decision-making process exhibit the
highest error rates, while comparatively lower error rates are observed at the
beginning and end phases of the process. Figure 4: The comparative analysis of
two key distributions: the frequency of operations in the training set (blue)
and the incidence of incorrect actions (red) encountered in the SSDNC dataset.
The y-axis represents error rates, while the x-axis enumerates various
operations. There is a notable correlation between these distributions,
indicating that the agent tends to favor operations it encountered more
frequently during training.

The trend of error rate throughout the expression generation process. To
account for expressions of different lengths, we normalize them by calculating
the percentage of completion. Notably, the highest error rates occur in
decisions made during the middle stage of the process, as depicted in Figure
3. A potential explanation for this is that initial decisions typically
involve binary operations, such as addition ($+$) and multiplication
($\times$). These are relatively straightforward, often following a general
strategy of expanding the tree at the outset. Conversely, the decisions made
towards the end of the process tend to be simpler, as the expression narrows
down and becomes more specific.

The distribution of incorrect selections. Our analysis revealed that the
frequency of errors closely reflects the distribution of operations in the
training set. Specifically, the agent is more inclined to select operations
that it encountered more frequently during training. This pattern is clearly
illustrated in Figure 4. This finding indicates a potential imbalance in the
design of the training set, a consideration that should be addressed in future
works.

Figure 5: Distribution of incorrect decision choices: These charts detail the
choices made during instances of incorrect decision-making. Each segment of
the pie chart corresponds to an incorrect choice related to a specific correct
operation. The size of each segment reflects the proportion of these incorrect
choices, with their total sum equaling $1$. The title of each chart identifies
the correct operation that was targeted in that specific instance.

Analysis of specific error types. We observed that the agent often struggles
with determining the exact values of constants. Although it correctly
identifies the necessity to include a constant, the chosen value frequently
turns out to be incorrect, as illustrated in Figure 5 e) to i). This issue
arises because the model does not apply coefficient fitting for constants;
instead, it treats each discrete value as a separate action. Additionally, the
agent commonly confuses sine and cosine functions. This confusion likely stems
from their similar mathematical properties, as shown in Figures 5 c) and d).

Another notable observation concerns the square root operation. The agent
often confuses it with the power operation, as illustrated in Figure 5 g).
This confusion is somewhat understandable, considering that the square root is
equivalent to raising a number to the power of $\frac{1}{2}$. This pattern may
suggest that the agent has some grasp of the impact of these operations on
mathematical relationships. A similar confusion occurs with the ’division’
$\div$ operation, which is often mistaken for the ’power’ operation since
’division’ is equivalent to ’power -1’, as shown in Figure 5 a).

##  6 Validation of Online Searching

A significant contribution of our study is the capability of our proposed
pretrained model to integrate effectively with online RL methods. We
demonstrate that our pretrained model can be seamlessly combined with the
online reinforcement learning method described in Section 3.7 to effectively
guide the search process. As a case study, we focused on $170$ challenging
instances within the SSDNC dataset, where the agent fails to correctly recover
the skeleton, even with beam search implementation. These instances present
substantial challenges for online exploration and learning. It is noteworthy
that we limit the agent to a maximum of $50$ explorations per search, a number
significantly lower than that in other online search methods like DSO
(Petersen et al., 2019), which often necessitate over $10,000$ searches. We
perform the online search for each instance individually. After the online
search process, we observe a significant performance improvement in $74.7\%$
of the cases. The initial average $R^{2}$ for the $170$ cases is $0.328$.
Remarkably, the average absolute improvement in $R^{2}$ for these unsolved
expressions is $0.342$, a substantial increase of $104.5\%$. The results are
summarised in Table 2.

Methods | $R^{2}$
---|---
Sym-Q | 0.328
Sym-Q + online search | 0.670
Table 2: The average $R^{2}$ after online search without beam search strategy.

Moreover, after the online search, the agent successfully uncovers some of
underlying skeletons, a significant achievement considering the pretrained
Sym-Q’s inability to detect them, even with a $128$-size beam search, while
the agent accomplishes this with only $50$ exploration. For example, in case:
$\sin(x_{1}^{2}+x_{1}x_{2}+x_{2})$, where the initial output is
$\sin(x_{2}+x_{1}^{3})$ rewarded $0$, the agent uncovers the underlying
expression by trying different combinations of $x_{1}$ and $x_{2}$ to replace
term inside the sine function.

##  7 Discussion and Conclusion

In conclusion, this study marks a significant milestone in the field of
symbolic regression. In this work, we innovatively reformulate the symbolic
regression as a sequential decision-making task and introduce the Sym-Q model.
Sym-Q is trained via our proprosed offline RL method and can be also
seamlessly combined with online RL approaches. We experimentally demonstrate
the effectiveness of this novel problem formulation. Additionally, Sym-Q has
the capability to directly derive potential expressions for many unseen
equations with high accuracy. Through extensive benchmark evaluations, Sym-Q
has been shown to outperform other state-of-the-art algorithms in terms of
skeleton recovery rate and to exhibit superior fitting accuracy across most
benchmarks. Notably, Sym-Q performs superior performance in recovering the
true underlying equations, making it particularly promising for discovering
unknown mathematical and physical problems in the future.

Most importantly, the open-source pretrained model we introduce could serve as
an effective starting point for more efficient and effective online search.
Additionally, we demonstrate for the first time how to integrate this
pretrained model with online RL and show its effectiveness in guiding the
search process.

A potential future research direction involves designing a more comprehensive
dataset that encompasses higher-dimensional expressions, including Ordinary
Differential Equations (ODE) and Partial Differential Equations (PDE), for
broader pretraining. Additionally, there is still room for optimizing the
current MDP formulation. Our current framework represents just the beginning
of a promising journey in this field.

## Potential Broader Impact

This paper presents work whose goal is to advance the field of Machine
Learning. There are many potential societal consequences of our work, none
which we feel must be specifically highlighted here.

## References

  * Augusto & Barbosa (2000) Augusto, D. A. and Barbosa, H. J.  Symbolic regression via genetic programming.  In _Proceedings. Vol. 1. Sixth Brazilian symposium on neural networks_ , pp. 173–178. IEEE, 2000.
  * Bendinelli et al. (2023) Bendinelli, T., Biggio, L., and Kamienny, P.-A.  Controllable neural symbolic regression.  _arXiv preprint arXiv:2304.10336_ , 2023.
  * Biggio et al. (2021) Biggio, L., Bendinelli, T., Neitz, A., Lucchi, A., and Parascandolo, G.  Neural symbolic regression that scales.  In _International Conference on Machine Learning_ , pp. 936–945. PMLR, 2021.
  * Bladek & Krawiec (2019) Bladek, I. and Krawiec, K.  Solving symbolic regression problems with formal constraints.  In _Proceedings of the Genetic and Evolutionary Computation Conference_ , pp. 977–984, 2019.
  * Devlin et al. (2018) Devlin, J., Chang, M.-W., Lee, K., and Toutanova, K.  Bert: Pre-training of deep bidirectional transformers for language understanding.  _arXiv preprint arXiv:1810.04805_ , 2018.
  * Dubčáková (2011) Dubčáková, R.  Eureqa: software review, 2011.
  * Fletcher (2000) Fletcher, R.  _Practical methods of optimization_.  John Wiley & Sons, 2000.
  * Forrest (1993) Forrest, S.  Genetic algorithms: principles of natural selection applied to computation.  _Science_ , 261(5123):872–878, 1993.
  * Glantz et al. (2001) Glantz, S. A., Slinker, B. K., and Neilands, T. B.  _Primer of applied regression & analysis of variance, ed_, volume 654.  McGraw-Hill, Inc., New York, 2001.
  * Gustafson et al. (2005) Gustafson, S., Burke, E. K., and Krasnogor, N.  On improving genetic programming for symbolic regression.  In _2005 IEEE Congress on Evolutionary Computation_ , volume 1, pp. 912–919. IEEE, 2005.
  * Hendrycks & Gimpel (2016) Hendrycks, D. and Gimpel, K.  Gaussian error linear units (gelus).  _arXiv preprint arXiv:1606.08415_ , 2016.
  * Hopcroft et al. (2006) Hopcroft, J. E., Motwani, R., and Ullman, J. D.  Automata theory, languages, and computation.  _International Edition_ , 24(2):171–183, 2006.
  * Janner et al. (2021) Janner, M., Li, Q., and Levine, S.  Offline reinforcement learning as one big sequence modeling problem.  _Advances in neural information processing systems_ , 34:1273–1286, 2021.
  * Kamienny et al. (2022) Kamienny, P.-A., d’Ascoli, S., Lample, G., and Charton, F.  End-to-end symbolic regression with transformers.  _Advances in Neural Information Processing Systems_ , 35:10269–10281, 2022.
  * Keijzer (2003) Keijzer, M.  Improving symbolic regression with interval arithmetic and linear scaling.  In _European Conference on Genetic Programming_ , pp. 70–82. Springer, 2003.
  * Kostrikov et al. (2021) Kostrikov, I., Nair, A., and Levine, S.  Offline reinforcement learning with implicit q-learning.  _arXiv preprint arXiv:2110.06169_ , 2021.
  * Koza (1994) Koza, J. R.  Genetic programming as a means for programming computers by natural selection.  _Statistics and computing_ , 4:87–112, 1994.
  * Krawiec & Pawlak (2013) Krawiec, K. and Pawlak, T.  Approximating geometric crossover by semantic backpropagation.  In _Proceedings of the 15th annual conference on Genetic and evolutionary computation_ , pp. 941–948, 2013.
  * Kumar et al. (2020) Kumar, A., Zhou, A., Tucker, G., and Levine, S.  Conservative q-learning for offline reinforcement learning.  _Advances in Neural Information Processing Systems_ , 33:1179–1191, 2020.
  * Lample & Charton (2019) Lample, G. and Charton, F.  Deep learning for symbolic mathematics.  _arXiv preprint arXiv:1912.01412_ , 2019.
  * Lee et al. (2019) Lee, J., Lee, Y., Kim, J., Kosiorek, A., Choi, S., and Teh, Y. W.  Set transformer: A framework for attention-based permutation-invariant neural networks.  In _International conference on machine learning_ , pp. 3744–3753. PMLR, 2019.
  * Levine et al. (2020) Levine, S., Kumar, A., Tucker, G., and Fu, J.  Offline reinforcement learning: Tutorial, review, and perspectives on open problems.  _arXiv preprint arXiv:2005.01643_ , 2020.
  * Li et al. (2022) Li, W., Li, W., Sun, L., Wu, M., Yu, L., Liu, J., Li, Y., and Tian, S.  Transformer-based model for symbolic regression via joint supervised learning.  In _The Eleventh International Conference on Learning Representations_ , 2022.
  * Ma et al. (2022) Ma, X., Qin, C., You, H., Ran, H., and Fu, Y.  Rethinking network design and local geometry in point cloud: A simple residual mlp framework.  _arXiv preprint arXiv:2202.07123_ , 2022.
  * Matchev et al. (2022) Matchev, K. T., Matcheva, K., and Roman, A.  Analytical modeling of exoplanet transit spectroscopy with dimensional analysis and symbolic regression.  _The Astrophysical Journal_ , 930(1):33, 2022.
  * Mnih et al. (2015) Mnih, V., Kavukcuoglu, K., Silver, D., Rusu, A. A., Veness, J., Bellemare, M. G., Graves, A., Riedmiller, M., Fidjeland, A. K., Ostrovski, G., et al.  Human-level control through deep reinforcement learning.  _nature_ , 518(7540):529–533, 2015.
  * Mundhenk et al. (2021) Mundhenk, T. N., Landajuela, M., Glatt, R., Santiago, C. P., Faissol, D. M., and Petersen, B. K.  Symbolic regression via neural-guided genetic programming population seeding.  _arXiv preprint arXiv:2111.00053_ , 2021.
  * Petersen et al. (2019) Petersen, B. K., Landajuela, M., Mundhenk, T. N., Santiago, C. P., Kim, S. K., and Kim, J. T.  Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients.  _arXiv preprint arXiv:1912.04871_ , 2019.
  * Riedmiller et al. (2018) Riedmiller, M., Hafner, R., Lampe, T., Neunert, M., Degrave, J., Wiele, T., Mnih, V., Heess, N., and Springenberg, J. T.  Learning by playing solving sparse reward tasks from scratch.  In _International conference on machine learning_ , pp. 4344–4353. PMLR, 2018.
  * Schmidt & Lipson (2009) Schmidt, M. and Lipson, H.  Distilling free-form natural laws from experimental data.  _science_ , 324(5923):81–85, 2009.
  * Silver et al. (2016) Silver, D., Huang, A., Maddison, C. J., Guez, A., Sifre, L., Van Den Driessche, G., Schrittwieser, J., Antonoglou, I., Panneershelvam, V., Lanctot, M., et al.  Mastering the game of go with deep neural networks and tree search.  _nature_ , 529(7587):484–489, 2016.
  * Silver et al. (2017) Silver, D., Schrittwieser, J., Simonyan, K., Antonoglou, I., Huang, A., Guez, A., Hubert, T., Baker, L., Lai, M., Bolton, A., et al.  Mastering the game of go without human knowledge.  _nature_ , 550(7676):354–359, 2017.
  * Tenachi et al. (2023) Tenachi, W., Ibata, R., and Diakogiannis, F. I.  Deep symbolic regression for physics guided by units constraints: toward the automated discovery of physical laws.  _arXiv preprint arXiv:2303.03192_ , 2023.
  * Tian et al. (2020) Tian, Y., Wang, Q., Huang, Z., Li, W., Dai, D., Yang, M., Wang, J., and Fink, O.  Off-policy reinforcement learning for efficient and effective gan architecture search.  In _Computer Vision–ECCV 2020: 16th European Conference, Glasgow, UK, August 23–28, 2020, Proceedings, Part VII 16_ , pp. 175–192. Springer, 2020.
  * Udrescu & Tegmark (2020) Udrescu, S.-M. and Tegmark, M.  Ai feynman: A physics-inspired method for symbolic regression.  _Science Advances_ , 6(16):eaay2631, 2020.
  * Uy et al. (2011) Uy, N. Q., Hoai, N. X., O’Neill, M., McKay, R. I., and Galván-López, E.  Semantically-based crossover in genetic programming: application to real-valued symbolic regression.  _Genetic Programming and Evolvable Machines_ , 12:91–119, 2011.
  * Valipour et al. (2021) Valipour, M., You, B., Panju, M., and Ghodsi, A.  Symbolicgpt: A generative transformer model for symbolic regression.  _arXiv preprint arXiv:2106.14131_ , 2021.
  * Vastl et al. (2022) Vastl, M., Kulhánek, J., Kubalík, J., Derner, E., and Babuška, R.  Symformer: End-to-end symbolic regression using transformer-based architecture.  _arXiv preprint arXiv:2205.15764_ , 2022.
  * Vaswani et al. (2017) Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, Ł., and Polosukhin, I.  Attention is all you need.  _Advances in neural information processing systems_ , 30, 2017.
  * Vecerik et al. (2017) Vecerik, M., Hester, T., Scholz, J., Wang, F., Pietquin, O., Piot, B., Heess, N., Rothörl, T., Lampe, T., and Riedmiller, M.  Leveraging demonstrations for deep reinforcement learning on robotics problems with sparse rewards.  _arXiv preprint arXiv:1707.08817_ , 2017.
  * Wadekar et al. (2022) Wadekar, D., Thiele, L., Villaescusa-Navarro, F., Hill, J. C., Cranmer, M., Spergel, D. N., Battaglia, N., Anglés-Alcázar, D., Hernquist, L., and Ho, S.  Augmenting astrophysical scaling relations with machine learning: application to reducing the sz flux-mass scatter.  _arXiv preprint arXiv:2201.01305_ , 2022.

##  Appendix A Pseudo-code for Deployment

We outline the Sym-Q framework in Figure 1. Here we provide the inference
pseudo-code in Algorithm 11.

Algorithm 1 Symbolic regression in the test time

Observe the tensor $T$ contains data points

Initialize the expression tree and the initial one-hot tree matrix $M$ filling
with all $0$

Initialize the beam search buffer $B$ with size $2*L$, and the output buffer
$OB$ with size $L$

while Output buffer not full or beam search buffer emptied do

Sample all the potential actions from $Q_{\theta}(T,M)$ and put into the beam
search buffer $B$

Save and update the log probability of the trajectories after appending the
new actions

Rank the trajectories in the beam search buffer $B$ according to the log
probability

Move the terminated trajectories to the output buffer $OB$

Keep the top $2*L$ trajectories in the beam search buffer $B$

end while

Performance BFGS for all the trajectories in the output buffer $OB$ and
provide the closed-form expressions

##  Appendix B Neural Networks and Hyperparameters

The Sym-Q consists of four main components: the points encoder, the tree
encoder, the fusion layer, and the Q-network. For the points encoder, we apply
the Set Transformer architecture (Devlin et al., 2018; Biggio et al., 2021),
using the original publicly available implementation, which has been proved
feasible for large-scale symbolic regression tasks (Biggio et al., 2021;
Bendinelli et al., 2023). Referring to the notation of the original paper, our
encoder is formed by multiple Induced Set Attention Blocks (ISABs) and one
final component performing Pooling by Multihead Attention (PMA), which can
reduce the computation burden associated with the self-attention operation.
PMA allows us to aggregate the output of the encoder into trainable abstract
features representing a compressed representation of the input equation. For
the tree encoder, we opted for a regular Transformer encoder (Vaswani et al.,
2017), using the default PyTorch implementation. In this research, we adopt a
fully connected layer of 4096 units as the fusion layer. Similarly, for the
Q-network, we use a fully connected MLP with four hidden layers of 4096 units,
predicting the Q-value. GELU activation function (Hendrycks & Gimpel, 2016) is
used in all the hidden layers in the fusion layer and Q-network. The other
hyper-parameters can be found in Table 3.

Table 3: Hyperparameters for our models. Points encoder
---
Number of ISABs | 5
Hidden dimension | 512
Number of heads | 8
Number of PMA features | 100
Number of inducing points | 100
Tree encoder
Number of layers | 5
Hidden dimension | 512
Number of heads | 8
Training
Batch size | 512
Learning rate | 0.0001
Scale weight (w/ CL) $\alpha$ | 0.2
Temperature parameter $\tau$ | 0.07
Max norm of gradients | 1.0
Inference
Beam size | 128
Re-start times of BFGS | 20

##  Appendix C Efficient Conservative Offline Q-learning

Conservative Q-learning (CQL) (Kumar et al., 2020) addresses the value over-
estimation issue in offline reinforcement learning by minimizing the values of
state-actions out of the data distribution and maximizing the values over the
data distribution. In this work, considering the sparse reward characteristic
of symbolic regression, that only the termination of each trajectory can be
rewarded and most of the sub-optimal trajectories should be rewarded
relatively low, we introduce a modified version of CQL can effectively learn
with the optimal offline data:

$\displaystyle J(\theta)\;=$
$\displaystyle\;\mathbb{E}_{(s,a^{i},s^{\prime},r)\sim
D}\underbrace{\frac{1}{2}[Q_{\theta}(s_{t},a^{i})-(r+\gamma\max_{a^{*}}\hat{Q}_{\hat{\theta}}(s_{t+1},a^{*}))]^{2}}_{\text{term
I: TD error}}+$ (7)
$\displaystyle\underbrace{\frac{1}{2}\mathbb{E}_{a_{j}|j\neq i,1\leq j\leq
m}[Q_{\theta}(s,a^{j})-0]^{2}}_{\text{term II: conservative Q loss }},$

where $\theta$ is the weight parameters of $Q$, $D$ is the training dataset,
$\gamma$ is the discounted factor, $a^{i}$ is the demonstration action at
state $s$, and ${a_{j}|j\neq i,1\leq j\leq m}$ is the actions not in the
demonstration, which can be consider as suboptimial actions. The
$\hat{Q}_{\hat{\theta}}$ is the target Q-network that has the same structure
as $Q$ but updated through exponentially moving average of the $Q$ weights
$\theta$ to stabilize the learning (Mnih et al., 2015). As mentioned above,
the suboptimal actions mostly can not lead to correct expression and should be
rewarded low. Then, the second term, the conservative Q loss, minimize the Q
values to $0$, the lower bound of this task, which is also similar to the idea
of Risk seeking policy gradient in DSO (Petersen et al., 2019).

Since the training data, all the demonstrations lead to a perfect recovery of
skeleton, $(r+\gamma\max_{a^{*}}\hat{Q}_{\hat{\theta}}(s_{t+1},a^{*}))$ is
equal to $1$ given discounted factor $\gamma=1$, which is feasible in this
case due to the Q value is bounded between $[0,1]$. Then we can rewrite the
objective function to :

$\displaystyle J(\theta)\;=$ $\displaystyle\;\mathbb{E}_{(s,a^{i})\sim
D}\underbrace{\frac{1}{2}[Q_{\theta}(s,a^{i})-1]^{2}}_{\text{term I: Q
loss}}+$ (8) $\displaystyle\underbrace{\frac{1}{2}\mathbb{E}_{a_{j}|j\neq
i,1\leq j\leq m}[Q_{\theta}(s,a^{j})-0]^{2}}_{\text{term II: conservative Q
loss }}.$

The modified objective function can simplify the neural networks without
maintaining the target Q-network $\hat{Q}$. Then, we can further simplify the
objective function, and make it very close to multi-class cross-entropy loss
in the sense that it emphasizes the correct prediction (or action in this
case) while penalizing the others. This simplification is particularly
suitable for symbolic regression tasks where the action space is discrete and
the optimal action at each state is distinctly defined. Building upon the
connection between the conservative Q-learning loss and the cross-entropy
loss, we derive the following reformulation:

Maximizing Selected Action’s Q Value: Simplify the first term
$\frac{1}{2}[Q_{\theta}(s,a^{i})-1]^{2}$ of Equation 8 by maximizing
$Q_{\theta}(s,a^{i})$:

$\max Q_{\theta}(s,a^{i}).$ (9)

Normalizing via Softmax Function: Transform Equation 9 to a probabilistic
interpretation by applying the softmax function:

$p(a^{i}|s)=\frac{e^{Q_{\theta}(s,a^{i})}}{\sum_{j=1}^{m}e^{Q_{\theta}(s,a^{j})}}.$
(10)

This inherently manages the Q values of non-selected actions in the second
term of Equation 8.

Minimizing Negative Log Probability: Maximize $p(a^{i}|s)$ by minimizing its
negative log:

$-\log(p(a^{i}|s_{t}))=-\log\left(\frac{e^{Q_{\theta}(s,a^{i})}}{\sum_{j=1}^{m}e^{Q_{\theta}(s,a^{j})}}\right).$
(11)

New Objective Function: Incorporate the above steps into a new objective
function:

$J(\theta)=-\mathbb{E}_{(s,a^{i})\sim\mathcal{D}}\left[\log\left(\frac{e^{Q_{\theta}(s,a^{i})}}{\sum_{j=1}^{m}e^{Q_{\theta}(s,a^{j})}}\right)\right].$
(12)

The final objective function represents the log probability of the optimal
action, mirroring the focus of cross-entropy on the correct class. The softmax
function in the denominator normalizes the Q values across all possible
actions, turning them into a probability distribution.

By adopting this modified objective function, our approach aligns more closely
with the principles of traditional supervised learning, potentially making it
more intuitive and easier to implement, especially in environments where the
correct action at each state is clearly defined, such as symbolic regression.

##  Appendix D Details for Dataset Generation and Action Space

The training set, and SSDNC benchmark are generated following the guidelines
provided by Li et al. (2022). We sample each non-leaf node following the
unnormalized weighted distribution shown in Table 4. The database
approximately contains $100,000$ unique expressions’ skeletons. Then the
training set is generated by re-sampling $50$ numerical constants for each
skeleton. The SSDNC benchmark includes approximately $100$ unique expression
skeletons and $10$ re-sampled numerical constants for each skeleton. In this
work, at every time step $t$, we consider $30$ different operations, which is
listed below in Table 5.

Operation | Mathematical meaning | Unnormalized probability
---|---|---
add | $+$ | 10
mul | $\times$ | 10
sub | $-$ | 5
div | $\div$ | 5
sqrt | $\sqrt{}$ | 4
exp | $\exp$ | 4
ln | $\ln$ | 4
sin | $\sin$ | 4
cos | $\cos$ | 4
tan | $\tan$ | 4
pow2 | $(\cdot)^{2}$ | 4
pow3 | $(\cdot)^{3}$ | 4
pow4 | $(\cdot)^{4}$ | 2
pow5 | $(\cdot)^{5}$ | 1
Table 4: Unnormalised probabilities of unary and binary operators used by the dataset generator. Operation | Operation id
---|---
$x_{1}$ | $0$
$x_{2}$ | $1$
$c$ | $2$
abs | $3$
$+$ | $4$
$*$ | $5$
$/$ | $6$
$\sqrt{}$ | $7$
$\exp$ | $8$
$\log$ | $9$
$**$ | $10$
$\sin$ | $11$
$\cos$ | $12$
$\tan$ | $13$
$\sin^{-1}$ | $14$
$\cos^{-1}$ | $15$
$\tan^{-1}$ | $16$
$\sinh$ | $17$
$\cosh$ | $18$
$\tanh$ | $19$
$\coth$ | $20$
$-3$ | $21$
$-2$ | $22$
$-1$ | $23$
$0$ | $24$
$1$ | $25$
$2$ | $26$
$3$ | $27$
$4$ | $28$
$5$ | $29$
Table 5: Action space and corresponding id.

##  Appendix E Fitting Accuracy on Additional Benchmarks

We evaluate our method and current state-of-the-art approaches on the widely
used public benchmarks, including Nguyen benchmark (Uy et al., 2011), the
Constant benchmark, the Keijzer benchmark (Keijzer, 2003), R rationals
(Krawiec & Pawlak, 2013), the AI-Feynman database (Udrescu & Tegmark, 2020),
as well as the SSDNC dataset (Petersen et al., 2019). Nguyen was the main
benchmark used in Petersen et al. (2019). In our study, we consider all the
expressions with up to two variables. From the results in Table 6, our method
outperforms all baseline methods in terms of weighted average $R^{2}$ on six
benchmarks. The weighted average $R^{2}$ is calculated based on the number of
expressions in each benchmark. In Table 7-11, we show the name of the
benchmark and corresponding expressions.

|  | Sym-Q | SymPointMLP | NeSymRes | SymbolicGPT | DSO | GP
---|---|---|---|---|---|---|---
Benchmark | Number of expressions | $R^{2}$ | $R^{2}$ | $R^{2}$ | $R^{2}$ | $R^{2}$ | $R^{2}$
Nguyen | 12 | 0.86482 | 0.99999 | 0.97538 | 0.64394 | 0.99489 | 0.89019
Constant | 8 | 0.86860 | 0.99998 | 0.84935 | 0.69433 | 0.99927 | 0.90842
Keijzer | 12 | 0.94392 | 0.98320 | 0.97500 | 0.59457 | 0.96928 | 0.90082
R | 3 | 0.99999 | 0.99999 | 0.84935 | 0.71093 | 0.97298 | 0.83198
AI-Feynman | 15 | 0.99975 | 0.99999 | 0.99999 | 0.64682 | 0.99999 | 0.92242
SSDNC | 963 | 0.95135 | 0.94782 | 0.85792 | 0.74585 | 0.93198 | 0.88913
Weighted avg. | 1013 | 0.95044 | 0.95020 | 0.86271 | 0.74087 | 0.93483 | 0.88976
Table 6: Results comparing our Sym-Q with state-of-the-art methods on several benchmarks. All methods use the beam search strategy with the beam size equaling 128. We report the average value of $R^{2}$ for each benchmark and the weighted average $R^{2}$ on six benchmarks. Name | Expression
---|---
Nguyen-1 | $x^{3}+x^{2}+x$
Nguyen-2 | $x^{4}+x^{3}+x^{2}+x$
Nguyen-3 | $x^{5}+x^{4}+x^{3}+x^{2}+x$
Nguyen-4 | $x^{6}+x^{5}+x^{4}+x^{3}+x^{2}+x$
Nguyen-5 | $\sin(x^{2})\cos(x)-1$
Nguyen-6 | $\sin(x)+\sin(x+x^{2})$
Nguyen-7 | $\ln(x)+\ln(x^{2}+1)$
Nguyen-8 | $\sqrt{x}$
Nguyen-9 | $\sin(x)+\sin(y)$
Nguyen-9 | $\sin(x)+\sin(y)$
Nguyen-10 | $2\sin(x)\cos(y)$
Nguyen-11 | $x^{y}$
Nguyen-12 | $x^{4}-x^{3}+\frac{1}{2}y^{2}-y$
Table 7: Benchmark functions in Nguyen. Name | Expression
---|---
Constant-1 | $3.39x^{3}+2.12x^{2}+1.78x$
Constant-2 | $\sin(x^{2})\cos(x)-0.75$
Constant-3 | $\sin(1.5x)\cos(0.5y)$
Constant-4 | $2.7x^{y}$
Constant-5 | $\sqrt{1.23x}$
Constant-6 | $x^{0.423}$
Constant-7 | $2\sin(1.3x)\cos(y)$
Constant-8 | $\ln(x+1.4)+\ln(x^{2}+1.3)$
Table 8: Benchmark functions in Constant. Name | Expression
---|---
Keijzer-3 | $0.3x\sin(2\pi x)$
Keijzer-4 | $x^{3}\exp(-x)\cos(x)\sin(x)(\sin(x^{2})\cos(x)-1))$
Keijzer-6 | $\frac{x(x+1)}{2}$
Keijzer-7 | $\ln{x}$
Keijzer-8 | $\sqrt{x}$
Keijzer-9 | $\ln{x+\sqrt{x^{2}+1}}$
Keijzer-10 | $x^{y}$
Keijzer-11 | $xy+\sin{(x-1)(y-1)}$
Keijzer-12 | $x^{4}-x^{3}+\frac{y^{2}}{2}-y$
Keijzer-13 | $6\sin{x}\cos{y}$
Keijzer-14 | $\frac{8}{2+x^{2}+y^{2}}$
Keijzer-15 | $\frac{x^{3}}{5}+\frac{y^{3}}{2}-y-x$
Table 9: Benchmark functions in Keijzer. Name | Expression
---|---
R-1 | $\frac{(x+1)^{3}}{x^{2}-x+1}$
R-2 | $\frac{x^{5}-3x^{3}+1}{x^{2}+1}$
R-3 | $\frac{x^{6}+x^{5}}{x^{4}+x^{3}+x^{2}+x}$
Table 10: Benchmark functions in R. Name | Expression
---|---
Feynman-1 | $\frac{\exp{\frac{-x^{2}}{2}}}{\sqrt{2\pi}}$
Feynman-2 | $\frac{\exp{\frac{-(xy^{-1})^{2}}{2}}}{\sqrt{2\pi}y}$
Feynman-3 | $xy$
Feynman-4 | $xy$
Feynman-5 | $\frac{1}{2}xy^{2}$
Feynman-6 | $\frac{x}{y}$
Feynman-7 | $\frac{\sin{x}}{\sin{y}}$
Feynman-8 | $\frac{x}{y}$
Feynman-9 | $\frac{xy}{2\pi}$
Feynman-10 | $1.5xy$
Feynman-11 | $\frac{x}{4\pi y^{2}}$
Feynman-12 | $\frac{xy^{2}}{x}$
Feynman-13 | $xy^{2}$
Feynman-14 | $\frac{x}{2(1+y)}$
Feynman-15 | $\frac{xy}{2\pi}$
Table 11: Benchmark functions in AI-Feynman.


11institutetext: Algorithms and Complexity Group, TU Wien, Vienna, Austria
11email<EMAIL_ADDRESS>
# Splitting Plane Graphs to Outerplanarity††thanks: Anaïs Villedieu is
supported by the Austrian Science Fund (FWF) under grant P31119.

Martin Gronemann Martin Nöllenburg Anaïs Villedieu

###### Abstract

Vertex splitting replaces a vertex by two copies and partitions its incident
edges amongst the copies. This problem has been studied as a graph editing
operation to achieve desired properties with as few splits as possible, most
often planarity, for which the problem is -hard. Here we study how to minimize
the number of splits to turn a plane graph into an outerplane one. We tackle
this problem by establishing a direct connection between splitting a plane
graph to outerplanarity, finding a connected face cover, and finding a
feedback vertex set in its dual. We prove -completeness for plane biconnected
graphs, while we show that a polynomial-time algorithm exists for maximal
planar graphs. Finally, we provide upper and lower bounds for certain families
of maximal planar graphs.

###### Keywords:

vertex splitting outerplanarity feedback vertex set.

##  1 Introduction

Graph editing problems are fundamental problems in graph theory. They define a
set of basic operations on a graph $G$ and ask for the minimum number of these
operations necessary in order to turn $G$ into a graph of a desired target
graph class $\mathcal{G}$ [34, 29, 42, 24]. For instance, in the Cluster
Editing problem [38] the operations are insertions or deletions of individual
edges and the target graph class are cluster graphs, i.e., unions of vertex-
disjoint cliques. In graph drawing, a particularly interesting graph class are
planar graphs, for which several related graph editing problems have been
studied, e.g., how many vertex deletions are needed to turn an arbitrary graph
into a planar one [32] or how many vertex splits are needed to obtain a planar
graph [23, 16]. In this paper, we are interested in the latter operation:
vertex splitting. A _vertex split_ creates two copies of a vertex $v$,
distributes its edges among these two copies and then deletes $v$ from $G$.

Further, we are translating the graph editing problem into a more geometric or
topological drawing editing problem. This means that we apply the splitting
operations not to the vertices of an abstract graph, but to the vertices of a
planar graph drawing, or more generally to a planar embedded (or _plane_)
graph. In a plane graph, each vertex has an induced cyclic order of incident
edges, which needs to be respected by any vertex split in the sense that we
must split its cyclic order into two contiguous intervals, one for each of the
two copies. From a different perspective, the two faces that serve as the
separators of these two edge intervals are actually merged into a single face
by the vertex split.

Finally, we consider outerplanar graphs as the target graph class. Thus, we
want to apply a minimum number of vertex splits to a plane graph $G$, which
merge a minimum number of faces in order to obtain an outerplanar embedded
graph $G^{\prime}$, where all vertices are incident to a single face, called
the _outer face_. We denote this minimum number of splits as the _outerplane
splitting number_ $\operatorname{osn}(G)$ of $G$ (see Fig. 1). Outerplanar
graphs are a prominent graph class in graph drawing (see, e.g., [7, 17, 27,
28]) as well as in graph theory and graph algorithms more generally (e.g.,
[10, 18, 31]). For instance, outerplanar graphs admit planar circular layouts
or 1-page book embeddings [5]. Additionally, outerplanar graphs often serve as
a simpler subclass of planar graphs with good algorithmic properties. For
instance, they have treewidth 2 and their generalizations to $k$-outerplanar
graphs still have bounded treewidth [9, 6], which allows for polynomial-time
algorithms for -complete problems that are tractable for such bounded-
treewidth graphs. This, in turn, can be used to obtain a PTAS for these
problems on planar graphs [4].

We are now ready to define our main computational problem as follows.

######  Problem 1 (Outerplane Splitting Number)

Given a plane biconnected graph $G=(V,E)$ and an integer $k$, can we transform
$G$ into an outerplane graph $G^{\prime}$ by applying at most $k$ vertex
splits to $G$?

(a)

(b)

Figure 1: (a) An instance of Outerplane Splitting Number, where the colored
vertices will be split; (b) resulting outerplane graph after the minimum 3
splits.

#### Contributions.

In this paper, we introduce the above problem Outerplane Splitting Number. We
start by showing the key property for our subsequent results, namely that
(minimum) sets of vertex splits to turn a plane biconnected graph $G$ into an
outerplane one correspond to (minimum) connected face covers in $G$ (Section
2), which in turn are equivalent to (minimum) feedback vertex sets in the dual
graph of $G$. Using this tool we then show that for general plane biconnected
graphs Outerplane Splitting Number is -complete (Section 3), whereas for
maximal planar graphs we can solve it in polynomial time (Section 4). Finally,
we provide upper and lower bounds on the outerplane splitting number for
maximal planar graphs (Section 5).

#### Related Work.

Splitting numbers have been studied mostly for abstract (non-planar) graphs
with the goal of turning them into planar graphs. The Planar Splitting Number
problem is -complete in general [16], but exact splitting numbers are known
for complete and complete bipartite graphs [21, 23], as well as for the 4-cube
[15]. For two-layer drawings of general bipartite graphs, the problem is still
-complete, but [2] when parametrized by the number of split vertices. It has
also been studied for other surfaces such as the torus [19] and the projective
plane [20]. Another related concept is the split thickness of a graph $G$ (or
its folded covering number [26]), which is the smallest $k$ such that $G$ can
be transformed into a planar graph by applying at most $k$ splits per vertex.
Recognizing graphs with split thickness $2$ is -hard, but there is a constant-
factor approximation algorithm and a fixed-parameter algorithm for graphs of
bounded treewidth [14]. Recently, the complexity of the embedded splitting
number problem of transforming non-planar graph drawings into plane ones has
been investigated [35]. Beyond the theoretical investigations of splitting
numbers and planarity, there are also applied work in graph drawing making use
of vertex splitting to untangle edges [41] or to improve layout quality for
community exploration [3, 22].

Regarding vertex splitting for achieving graph properties other than
planarity, Trotter and Harary [39] studied vertex splitting to turn a graph
into an interval graph. Paik et al. [36] considered vertex splitting to remove
long paths in directed acyclic graphs and Abu-Khzam et al. [1] studied
heuristics using vertex splitting for a cluster editing problem.

(a)

(b)

(c)

(d)

Figure 2: (a) Two touching faces $f_{1},f_{2}$ with a common vertex $v$ on
their boundary. (b) Result of the split of $v$ with respect to $f_{1},f_{2}$
joining them into a new face $f$. (c-d) Merging 4 faces $f_{1},\ldots,f_{4}$
covering a single vertex $v$ with 3 splits.

#### Preliminaries.

The key concept of our approach is to merge a set of faces of a given plane
graph $G=(V,E)$ with vertex set $V=V(G)$ and edge set $E=E(G)$ into one big
face which is incident to all vertices of $G$. Hence, the result is
outerplanar. The idea is that if two faces $f_{1}$ and $f_{2}$ share a vertex
$v$ on their boundary (we say $f_{1}$ and $f_{2}$ _touch_ , see Fig. 2(a)),
then we can split $v$ into two new vertices $v_{1},v_{2}$. In this way, we are
able to create a narrow gap, which merges $f_{1},f_{2}$ into a bigger face $f$
(see Fig. 2(b)). With this in mind, we formally define an _embedding-
preserving split_ of a vertex $v$ w.r.t. two incident faces $f_{1}$ and
$f_{2}$. We construct a new plane graph $G^{\prime}=(V^{\prime},E^{\prime})$
with $V^{\prime}=V\setminus\\{v\\}\cup\\{v_{1},v_{2}\\}$. Consider the two
neighbors of $v$ both incident to $f_{1}$ and let $w_{1}$ be the second
neighbor in clockwise order. Similarly, let $w_{i}$ be the second vertex
adjacent to $v$ and incident to $f_{2}$. We call $w_{d}$ the vertex preceding
$w_{1}$ in the cyclic ordering or the neighbors, with $d$ being the degree of
$v$, see Fig. 2(a). Note that while $w_{1}=w_{i-1}$ and $w_{i}=w_{d}$ is
possible, $w_{d}\neq w_{1}$ and $w_{i-1}\neq w_{i}$. For the set of edges, we
now set
$E^{\prime}=E\setminus\\{(v,w_{1}),\ldots,(v,w_{d})\\}\cup\\{(v_{2},w_{1}),\ldots,(v_{2},w_{i-1})\\}\cup\\{(v_{1},w_{i}),\ldots,(v_{1},w_{d})\\}$
and assume that they inherit their embedding from $G$. From now on we refer to
this operation simply as a _split_ or when $f_{1},f_{2}$ are clear from the
context, we may refer to _merging_ the two faces at $v$. The vertices
$v_{1},v_{2}$ introduced in place of $v$ are called _copies_ of $v$. If a copy
$v_{i}$ of a vertex $v$ is split again, then any copy of $v_{i}$ is also
called a copy of the original vertex $v$.

We can now reformulate the task of using as few splits as possible. Our
objective is to find a set of faces $S$ that satisfies two conditions. (1)
Every vertex in $G$ has to be on the boundary of at least one face $f\in S$,
that is, the faces in $S$ _cover_ all vertices in $V$.111Testing whether such
$S$ with $|S|\leq k$ exists, is the -complete problem Face Cover [8]. And (2)
for every two faces $f,f^{\prime}\in S$ there exists a set of faces
$\\{f_{1},\dots,f_{k}\\}\subseteq S$ such that
$f=f_{1},\ldots,f_{k}=f^{\prime}$, and $f_{i}$ touches $f_{i+1}$ for $1\leq
i<k$. In other words, $S$ is connected in terms of touching faces. We now
introduce the main tool in our constructions that formalizes this concept.

##  2 Face-Vertex Incidence Graph

Let $G=(V,E)$ be a plane biconnected graph and $F$ its set of faces. The
_face-vertex incidence graph_ is defined as $H=(V\cup F,E_{H})$ and contains
the edges $E_{H}=\\{(v,f)\in V\times F:v\text{ is on the boundary of }f\\}$.
Graph $H$ is by construction bipartite and we assume that it is plane by
placing each vertex $f\in F$ into its corresponding face in $G$.

###### Definition 1

Let $G$ be a plane biconnected graph, let $F$ be the set of faces of $G$, and
let $H$ be its face-vertex incidence graph. A _face cover_ of $G$ is a set
$S\subseteq F$ of faces such that every vertex $v\in V$ is incident to at
least one face in $S$. A face cover $S$ of $G$ is a _connected_ face cover if
the induced subgraph $H[S\cup V]$ of $S\cup V$ in $H$ is connected.

We point out that the problem of finding a connected face cover is not
equivalent to the Connected Face Hitting Set Problem [37], where a connected
set of vertices incident to every face is computed. We continue with two
lemmas that are concerned with merging multiple faces at the same vertex (Fig.
2(c)).

###### Lemma 1

Let $G$ be a plane biconnected graph and $S\subseteq F$ a subset of the faces
$F$ of $G$ that all have the vertex $v\in V$ on their boundary. Then $|S|-1$
splits are sufficient to merge the faces of $S$ into one.

###### Proof

Let $f_{1},\ldots,f_{k}$ with $k=|S|$ be the faces of $S$ in the clockwise
order as they appear around $v$ ($f_{1}$ chosen arbitrarily). We iteratively
merge $f_{1}$ with $f_{i}$ for $2\leq i\leq k$, which requires in total
$|S|-1$ splits (see Fig. 2(c) and Fig. 2(d)). ∎

###### Lemma 2

Let $G$ be a plane biconnected graph and let $S$ be a connected face cover of
$G$. Then $|S|-1$ splits are sufficient to merge the faces of $S$ into one.

###### Proof

Let $H^{\prime}=H[S\cup V]$ and compute a spanning tree $T$ in $H^{\prime}$.
For every vertex $v\in V(T)\cap V(G)$, we apply Lemma 1 with the face set
$F^{\prime}(v)=\\{f\in S\cap V(T)\;|\;(v,f)\in E(T)\\}$. We root the tree at
an arbitrary face $f^{\prime}\in S$, which provides a hierarchy on the
vertices and faces in $T$. Every vertex $v\in V(T)\cap V(G)$ requires by Lemma
1 $|F^{\prime}(v)|-1$ splits. Note that that for all leaf vertices in $T$,
$|F^{\prime}(v)|=1$, i.e., they will not be split. Each split is charged to
the children of $v$ in $T$. Since $H$ is bipartite, so is $T$. It follows that
every face $f\in S\setminus\\{f^{\prime}\\}$ is charged exactly once by its
parent, thus $|S|-1$ splits suffice. ∎

###### Lemma 3

Let $G$ be a plane biconnected graph and $\sigma$ a sequence of $k$ splits to
make $G$ outerplane. Then $G$ has a connected face cover of size $k+1$.

###### Proof

Since by definition applying $\sigma$ to $G$ creates a single big face that is
incident to all vertices in $V(G)$ by iteratively merging pairs of original
faces defining a set $S\subseteq F$, it is clear that $S$ is a face cover of
$G$ and since the result of the vertex splits and face merges creates a single
face, set $S$ must also be connected. ∎

As a consequence of Lemmas 2 and 3 we obtain that Outerplane Splitting Number
and computing a minimum connected face cover are equivalent.

###### Theorem 2.1

Let $G$ be a plane biconnected graph. Then $G$ has outerplane splitting number
$k$ if and only if it has a connected face cover of size $k+1$.

##  3 -completeness

In this section, we prove that finding a connected face cover of size $k$ (and
thus Outerplane Splitting Number) is -complete. The idea is to take the dual
of a planar biconnected Vertex Cover instance and subdivide every edge once
(we call this an _all-1-subdivision_). Note that the all-1-subdivision of a
graph $G$ corresponds to its vertex-edge incidence graph and the
all-1-subdivision of the dual of $G$ corresponds to the face-edge incidence
graph of $G$. A connected face cover then corresponds to a vertex cover in the
original graph, and vice versa. The following property greatly simplifies the
arguments regarding Definition 1.

###### Property 1

Let $G^{\prime}$ be an all-1-subdivision of a biconnected planar graph $G$ and
$S$ a set of faces that cover $V(G^{\prime})$. Then $S$ is a connected face
cover of $G^{\prime}$.

###### Proof

Let $H$ be the all-1-subdivision of the dual of $G$, and assume to the
contrary that the induced subgraph $H^{\prime}=H[S\cup V(G)]$ is not
connected. Then there exists an edge $(u,v)\in E(G)$ such that $u$ and $v$ are
in different connected components in $H^{\prime}$. Let $w$ be the subdivision
vertex of $(u,v)$ in $G^{\prime}$. As a subdivision vertex, $w$ is incident to
only two faces, one of which, say $f$, must be contained in $S$. But $f$ is
also incident to $u$ and $v$ and hence $u$ and $v$ are in the same component
of $H^{\prime}$ via face $f$, a contradiction. Hence $H^{\prime}$ is connected
and $S$ is a connected face cover of $G^{\prime}$. ∎

(a) primal (black) and
dual graph (gray)

(b) subdivision of
the dual graph

(c) vertex cover (red) in
the primal graph

Figure 3: Link between the primal graph $G$, its vertex cover, the dual $D$
and its subdivision $D^{*}$.

The proof of the next theorem is very similar to the reduction of Bienstock
and Monma to show -completeness of Face Cover [8]; due to differences in the
problem definitions, such as the connectivity of the face cover and whether
the input graph is plane or not, we provide the full reduction for the sake of
completeness.

###### Theorem 3.1

Deciding whether a plane biconnected graph $G$ has a connected face cover of
size at most $k$ is -complete.

###### Proof

Clearly the problem is in $\NP$. To prove hardness, we first introduce some
notation. Let $G$ be a plane biconnected graph and $D$ the corresponding dual
graph. Furthermore, let $D^{*}$ be the all-1-subdivision of $D$. We prove now
that a connected face cover $S^{*}$ of size $k$ in $D^{*}$ is in a one-to-one
correspondence with a vertex cover $S$ of size $k$ in $G$ (see Fig. 3). More
specifically, we show that the dual vertices of the faces of $S^{*}$ that form
a connected face cover in $D^{*}$, are a vertex cover for $G$ and vice versa.
The reduction is from the -complete Vertex Cover problem in biconnected planar
graphs in which all vertices have degree 3 (cubic graphs) [33].

####  Connected Face Cover $\Rightarrow$ Vertex Cover:

Let $G$ be such a biconnected plane Vertex Cover instance. Assume we have a
connected face cover $S^{*}$ with $|S^{*}|=k$ for $D^{*}$. Note that the faces
of $D^{*}$ correspond to the vertices in $G$. We claim that the faces $S^{*}$,
when mapped to the corresponding vertices $S\subseteq V(G)$ are a vertex cover
for $G$. Assume otherwise, that is, there exists an edge $e^{*}\in E(G)$ that
has no endpoints in $S$. However, $e^{*}$ has a dual edge $e\in E(D)$ and
therefore a subdivision vertex $v_{e}\in V(D^{*})$. Hence, there is a face
$f\in S^{*}$ that has $v_{e}$ on its boundary by definition of connected face
cover. And when mapped to $D$, $f$ has $e$ on its boundary, which implies that
the primal edge $e^{*}$ has at least one endpoint in $S^{*}$; a contradiction.

####  Vertex Cover $\Rightarrow$ Connected Face Cover:

To prove that a vertex cover $S$ induces a connected face cover $S^{*}$ in
$D^{*}$, we have to prove that $S^{*}$ covers all vertices and the induced
subgraph in the face-vertex incidence graph $H$ is connected. We proceed as in
the other direction. $S$ covers all edges in $E(G)$, thus every edge $e\in
E(D)$ is bounded by at least one face of $S^{*}$. Hence, every subdivision
vertex in $V(D^{*})$ is covered by a face of $S^{*}$. Furthermore, every
vertex in $D^{*}$ is adjacent to a subdivision vertex, thus, also covered by a
face in $S^{*}$. Since $S^{*}$ is covering all vertices, we obtain from
Property 1 that $S^{*}$ is a connected face cover. ∎

##  4 Feedback Vertex Set Approach

A _feedback vertex set_ $S^{\circ}\subset V(G)$ of a graph $G$ is a vertex
subset such that the induced subgraph $G[V(G)\setminus S^{\circ}]$ is acyclic.
We show here that finding a connected face cover $S$ of size $k$ for a plane
biconnected graph $G$ is equivalent to finding a feedback vertex set
$S^{\circ}\subset V(D)$ of size $k$ in the dual graph $D$ of $G$. The _weak
dual_ , i.e., the dual without a vertex for the outer face, of an outerplanar
graph is a forest. Thus we must find the smallest number of splits in $G$
which transform $D$ into a forest. In other words, we must must break all the
cycles in $D$, and hence all of the vertices in the feedback vertex set
$S^{\circ}$ of $D$ correspond to the faces of $G$ that should be merged
together (see Fig. 4).

(a) primal (gray) and
dual graph (black)

(b) spanning tree (black) of
the connected face cover

(c) outerplanar graph after
splitting along the tree

Figure 4: The connected face cover (blue) is a feedback vertex set (red) in
the dual.

###### Property 2

Let $H$ be the face-vertex incidence graph of a plane biconnected graph $G$
and let $S^{\circ}$ be a feedback vertex set in the dual $D$ of $G$. Then
$S^{\circ}$ induces a connected face cover $S$ in $G$.

###### Proof

We need to show that $S^{\circ}$ is a face cover and that it is connected.
First, assume there is a vertex $v\in V(G)$ of degree $\deg(v)=d$ that is not
incident to a vertex in $S^{\circ}$, i.e., a face of $G$. Since $G$ is
biconnected, $v$ is incident to $d$ faces $f_{1},\dots,f_{d}$, none of which
is contained in $S^{\circ}$. But then $D[V(D)\setminus S^{\circ}]$ has a cycle
$(f_{1},\dots,f_{d})$, a contradiction.

Next, we define $\overline{S^{\circ}}=V(D)\setminus S^{\circ}$ as the
complement of the feedback vertex set $S^{\circ}$ in $D$. Assume that $H[V\cup
S^{\circ}]$ has at least two separate connected components $C_{1},C_{2}$. Then
there must exist a closed curve in the plane separating $C_{1}$ from $C_{2}$,
which avoids the faces in $S^{\circ}$ and instead passes through a sequence
$(f_{1},\dots,f_{\ell})$ of faces in $\overline{S^{\circ}}$, where each pair
$(f_{i},f_{i+1})$ for $i\in\\{1,\dots,\ell-1\\}$ as well as $(f_{\ell},f_{1})$
are adjacent in the dual $D$. Again this implies that there is a cycle in
$D[V(D)\setminus S^{\circ}]$, a contradiction. Thus $S^{\circ}$ is a connected
face cover. ∎

###### Theorem 4.1

A plane biconnected graph $G$ has outerplane splitting number $k$ if and only
if its dual $D$ has a minimum feedback vertex set of size $k+1$.

###### Proof

Let $S^{\circ}$ be a minimum feedback vertex set of the dual $D$ of $G$ with
cardinality $|S^{\circ}|=k+1$ and let $H$ be the face-vertex incidence graph
of $G$. We know from Property 2 that $H^{\prime}=H[V(G)\cup S^{\circ}]$ is
connected and hence $S^{\circ}$ induces a connected face cover $S$ with
$|S|=k+1$. Then by Lemma 2 $G$ has $\operatorname{osn}(G)\leq k$.

Let conversely $\sigma$ be a sequence of $k$ vertex splits that turn $G$ into
an outerplane graph $G^{\prime}$ and let $F$ be the set of faces of $G$. By
Lemma 3 we obtain a connected face cover $S$ of size $k+1$ consisting of all
faces that are merged by $\sigma$. The complement $\overline{S}=F\setminus S$
consists of all faces of $G$ that are not merged by the splits in $\sigma$ and
thus are the remaining (inner) faces of the outerplane graph $G^{\prime}$.
Since $G^{\prime}$ is outerplane and biconnected, $\overline{S}$ is the vertex
set of the weak dual of $G^{\prime}$, which must be a tree. Hence $S$ is a
feedback vertex set in $D$ of size $k+1$ and the minimum feedback vertex set
in $D$ has size at most $k+1$. ∎

Since all faces in a maximal planar graph are triangles, the maximum vertex
degree of its dual is 3. Thus, we can apply the polynomial-time algorithm of
Ueno et al. [40] to this dual, which computes the minimum feedback vertex set
in graphs of maximum degree 3 by reducing the instance to polynomial-solvable
matroid parity problem instance, and obtain

###### Corollary 1

We can solve Outerplane Splitting Number for maximal planar graphs in
polynomial time.

Many other existing results for feedback vertex set extend to Outerplane
Splitting Number, e.g., it has a kernel of size $13k$ [11] and admits a [13].

##  5 Lower and Upper Bounds

In this section we provide some upper and lower bounds on the outerplane
splitting number in certain maximal planar graphs.

###  5.1 Upper Bounds

Based on the equivalence of Theorem 4.1 we obtain upper bounds on the
outerplane splitting number from suitable upper bounds on the feedback vertex
set problem, which has been studied for many graph classes, among them cubic
graphs [12]. Liu and Zhao [30] showed that cubic graphs $G=(V,E)$ of girth at
least four (resp., three) have a minimum feedback vertex set of size at most
$\frac{|V|}{3}$ (resp., $\frac{3|V|}{8}$). Kelly and Liu [25] showed that
connected planar subcubic graphs of girth at least five have a minimum
feedback vertex set of size at most $\frac{2|V|+2}{7}$. Recall that the girth
of a graph is the length of its shortest cycle.

###### Proposition 1

The outerplane splitting number of a maximal planar graph $G=(V,E)$ of minimum
degree (i) 3, (ii) 4, and (iii) 5, respectively, is at most (i)
$\frac{3|V|-10}{4}$, (ii) $\frac{2|V|-7}{3}$, and (iii) $\frac{4|V|-13}{7}$,
respectively.

###### Proof

Maximal planar graphs with $n=|V|$ vertices have $2n-4$ faces. So the
corresponding dual graphs have $2n-4$ vertices. Moreover, since the degree of
a vertex in $G$ corresponds to the length of a facial cycle in the dual,
graphs with minimum vertex degree 3, 4, or 5 have duals with girth 3, 4, or 5,
respectively. So if the minimum degree in $G$ is 3, we obtain an upper bound
on the feedback vertex set of $(3n-6)/4$; if the minimum degree is 4, the
bound is $(2n-4)/3$; and if the minimum degree is 5, the bound is $(4n-6)/7$.
The claim then follows from Theorem 4.1. ∎

###  5.2 Lower Bounds

We first provide a generic lower bound for the outerplane splitting number of
maximal planar graphs. Let $G$ be an $n$-vertex maximal planar graph with
$2n-4$ faces. Each face is a triangle incident to three vertices. In a
minimum-size connected face cover $S^{*}$, the first face covers three
vertices. Due to the connectivity requirement, all other faces can add at most
two newly covered vertices. Hence we need at least $\frac{n-1}{2}$ faces in
any connected face cover. By Theorem 2.1 this implies that
$\operatorname{osn}(G)\geq\frac{n-3}{2}$.

###### Proposition 2

Any maximal planar graph $G$ has outerplane splitting number at least
$\frac{|V(G)|-3}{2}$.

Next, towards a better bound, we define a family of maximal planar graphs
$T_{d}=(V_{d},E_{d})$ of girth 3 for $d\geq 0$ that have outerplane splitting
number at least $\frac{2|V_{d}|-8}{3}$. The family are the complete planar
3-trees of depth $d$, which are defined recursively as follows. The graph
$T_{0}$ is the 4-clique $K_{4}$. To obtain $T_{d}$ from $T_{d-1}$ for $d\geq
1$ we subdivide each inner triangular face of $T_{d-1}$ into three triangles
by inserting a new vertex and connecting it to the three vertices on the
boundary of the face.

###### Proposition 3

The complete planar 3-tree $T_{d}$ of depth $d$ has outerplane splitting
number at least $\frac{2|V_{d}|-8}{3}$.

###### Proof

Each $T_{d}$ is a maximal planar graph with
$n_{d}=3+\sum_{i=0}^{d}3^{i}=\frac{3^{d+1}+5}{2}$ vertices. All $3^{d}$ leaf-
level vertices added into the triangular faces of $T_{d-1}$ in the last step
of the construction have degree 3 and are incident to three exclusive faces,
i.e., there is no face that covers more than one of these leaf-level vertices.
This immediately implies that any face cover of $T_{d}$, connected or not, has
size at least $3^{d}$. From $n_{d}=\frac{3^{d+1}+5}{2}$ we obtain
$d=\log_{3}\frac{2n_{d}-5}{3}$ and $3^{d}=\frac{2n_{d}-5}{3}$. Theorem 2.1
then implies that $\operatorname{osn}(T_{d})\geq\frac{2n_{d}-8}{3}$. ∎

##  6 Open Problems

We have introduced the Outerplane Splitting Number problem and established its
complexity for plane biconnected graphs. The most important open question
revolves around the embedding requirement. Splitting operations can be defined
more loosely and allow for any new embedding and neighborhood of the split
vertices. In general, it is also of interest to understand how the problem
differs when the input graph does not have an embedding at all, as in the
original splitting number problem. Since Outerplane Splitting Number can be
solved in polynomial time for maximal planar graphs but is hard for plane
biconnected graphs, there is a complexity gap to be closed when faces of
degree more than three are involved. Vertex splitting in graph drawings has so
far been studied to achieve planarity and outerplanarity. A natural extension
is to study it for other graph classes or graph properties.

## References

  * [1] Abu-Khzam, F.N., Barr, J.R., Fakhereldine, A., Shaw, P.: A greedy heuristic for cluster editing with vertex splitting. In: Proc. 4th International Conference on Artificial Intelligence for Industries (AI4I). pp. 38–41. IEEE (2021). https://doi.org/10.1109/AI4I51902.2021.00017
  * [2] Ahmed, R., Kobourov, S.G., Kryven, M.: An FPT algorithm for bipartite vertex splitting. In: Angelini, P., von Hanxleden, R. (eds.) Proc. 30th International Symposium on Graph Drawing and Network Visualization (GD’22). LNCS, Springer (2022). https://doi.org/10.48550/arXiv.2208.12898
  * [3] Angori, L., Didimo, W., Montecchiani, F., Pagliuca, D., Tappini, A.: Hybrid graph visualizations with chordlink: Algorithms, experiments, and applications. IEEE Trans. Vis. and Comput. Graph. 28(2), 1288–1300 (2022). https://doi.org/10.1109/TVCG.2020.3016055
  * [4] Baker, B.S.: Approximation algorithms for NP-complete problems on planar graphs. J. ACM 41(1), 153–180 (1994). https://doi.org/10.1145/174644.174650
  * [5] Bernhart, F., Kainen, P.C.: The book thickness of a graph. J. Comb. Theory, Ser. B 27(3), 320–331 (1979). https://doi.org/10.1016/0095-8956(79)90021-2
  * [6] Biedl, T.: On triangulating k-outerplanar graphs. Discret. Appl. Math. 181, 275–279 (2015). https://doi.org/10.1016/j.dam.2014.10.017
  * [7] Biedl, T.C.: Small drawings of outerplanar graphs, series-parallel graphs, and other planar graphs. Discret. Comput. Geom. 45(1), 141–160 (2011). https://doi.org/10.1007/s00454-010-9310-z
  * [8] Bienstock, D., Monma, C.L.: On the complexity of covering vertices by faces in a planar graph. SIAM J. Comput. 17(1), 53–76 (1988). https://doi.org/10.1137/0217004
  * [9] Bodlaender, H.L.: A partial $k$-arboretum of graphs with bounded treewidth. Theor. Comput. Sci. 209(1-2), 1–45 (1998). https://doi.org/10.1016/S0304-3975(97)00228-4
  * [10] Bodlaender, H.L., Fomin, F.V.: Approximation of pathwidth of outerplanar graphs. J. Algorithms 43(2), 190–200 (2002). https://doi.org/10.1016/S0196-6774(02)00001-9
  * [11] Bonamy, M., Kowalik, L.: A 13k-kernel for planar feedback vertex set via region decomposition. Theor. Comput. Sci. 645, 25–40 (2016). https://doi.org/10.1016/j.tcs.2016.05.031
  * [12] Bondy, J.A., Hopkins, G., Staton, W.: Lower bounds for induced forests in cubic graphs. Can. Math. Bull. 30(2), 193–199 (1987). https://doi.org/10.4153/CMB-1987-028-5
  * [13] Demaine, E., Hajiaghayi, M.: Bidimensionality: New connections between FPT algorithms and PTASs. In: Proc. 37th Annual ACM Symposium on Theory of Computing (STOC). pp. 590–601 (01 2005). https://doi.org/10.1145/1070432.1070514
  * [14] Eppstein, D., Kindermann, P., Kobourov, S.G., Liotta, G., Lubiw, A., Maignan, A., Mondal, D., Vosoughpour, H., Whitesides, S., Wismath, S.K.: On the planar split thickness of graphs. Algorithmica 80(3), 977–994 (2018). https://doi.org/10.1007/s00453-017-0328-y
  * [15] Faria, L., de Figueiredo, C.M.H., de Mendonça N., C.F.X.: The splitting number of the 4-cube. In: Lucchesi, C.L., Moura, A.V. (eds.) Proc. 3rd Latin American Symposium on Theoretical Informatics (LATIN). LNCS, vol. 1380, pp. 141–150. Springer (1998). https://doi.org/10.1007/BFb0054317
  * [16] Faria, L., de Figueiredo, C.M.H., de Mendonça N., C.F.X.: Splitting number is NP-complete. Discret. Appl. Math. 108(1), 65–83 (2001). https://doi.org/10.1016/S0166-218X(00)00220-1
  * [17] Frati, F.: Planar rectilinear drawings of outerplanar graphs in linear time. Computat. Geom. 103, 101854 (2022). https://doi.org/10.1016/j.comgeo.2021.101854
  * [18] Frederickson, G.N.: Searching among intervals and compact routing tables. Algorithmica 15(5), 448–466 (1996). https://doi.org/10.1007/BF01955044
  * [19] Hartsfield, N.: The toroidal splitting number of the complete graph k${}_{\mbox{n}}$. Discret. Math. 62(1), 35–47 (1986). https://doi.org/10.1016/0012-365X(86)90039-7
  * [20] Hartsfield, N.: The splitting number of the complete graph in the projective plane. Graphs Comb. 3(1), 349–356 (1987). https://doi.org/10.1007/BF01788557
  * [21] Hartsfield, N., Jackson, B., Ringel, G.: The splitting number of the complete graph. Graphs Comb. 1(1), 311–329 (1985). https://doi.org/10.1007/BF02582960
  * [22] Henry, N., Bezerianos, A., Fekete, J.: Improving the readability of clustered social networks using node duplication. IEEE Trans. Vis. Comput. Graph. 14(6), 1317–1324 (2008). https://doi.org/10.1109/TVCG.2008.141
  * [23] Jackson, B., Ringel, G.: The splitting number of complete bipartite graphs. Arch. Math. 42(2), 178–184 (1984). https://doi.org/10.1007/BF01772941
  * [24] Kant, G.: Augmenting outerplanar graphs. J. Algorithms 21(1), 1–25 (1996). https://doi.org/10.1006/jagm.1996.0034
  * [25] Kelly, T., Liu, C.: Minimum size of feedback vertex sets of planar graphs of girth at least five. Eur. J. Comb. 61, 138–150 (2017). https://doi.org/10.1016/j.ejc.2016.10.009
  * [26] Knauer, K.B., Ueckerdt, T.: Three ways to cover a graph. Discret. Math. 339(2), 745–758 (2016). https://doi.org/10.1016/j.disc.2015.10.023
  * [27] Lazard, S., Lenhart, W.J., Liotta, G.: On the edge-length ratio of outerplanar graphs. Theor. Comput. Sci. 770, 88–94 (2019). https://doi.org/10.1016/j.tcs.2018.10.002
  * [28] Lenhart, W., Liotta, G.: Proximity drawings of outerplanar graphs. In: North, S.C. (ed.) Proc. 4th International Symposium on Graph Drawing (GD). LNCS, vol. 1190, pp. 286–302. Springer (1996). https://doi.org/10.1007/3-540-62495-3_55
  * [29] Lewis, J.M., Yannakakis, M.: The node-deletion problem for hereditary properties is NP-complete. J. Comput. Syst. Sci. 20(2), 219–230 (1980). https://doi.org/10.1016/0022-0000(80)90060-4
  * [30] Liu, J., Zhao, C.: A new bound on the feedback vertex sets in cubic graphs. Discret. Math. 148(1-3), 119–131 (1996). https://doi.org/10.1016/0012-365X(94)00268-N
  * [31] Maheshwari, A., Zeh, N.: External memory algorithms for outerplanar graphs. In: Aggarwal, A., Rangan, C.P. (eds.) Proc. 10th International Symposium on Algorithms and Computation (ISAAC). LNCS, vol. 1741, pp. 307–316. Springer (1999). https://doi.org/10.1007/3-540-46632-0_31
  * [32] Marx, D., Schlotter, I.: Obtaining a planar graph by vertex deletion. Algorithmica 62(3-4), 807–822 (2012). https://doi.org/10.1007/s00453-010-9484-z
  * [33] Mohar, B.: Face covers and the genus problem for apex graphs. J. Comb. Theory, Ser. B 82(1), 102–117 (2001). https://doi.org/10.1006/jctb.2000.2026
  * [34] Natanzon, A., Shamir, R., Sharan, R.: Complexity classification of some edge modification problems. Discret. Appl. Math. 113(1), 109–128 (2001). https://doi.org/10.1016/S0166-218X(00)00391-7
  * [35] Nöllenburg, M., Sorge, M., Terziadis, S., Villedieu, A., Wu, H., Wulms, J.: Planarizing graphs and their drawings by vertex splitting. In: Angelini, P., von Hanxleden, R. (eds.) Proc. 30th International Symposium on Graph Drawing and Network Visualization (GD’22). LNCS, Springer (2022), https://arxiv.org/abs/2202.12293
  * [36] Paik, D., Reddy, S.M., Sahni, S.: Vertex splitting in dags and applications to partial scan designs and lossy circuits. Int. J. Found. Comput. Sci. 9(4), 377–398 (1998). https://doi.org/10.1142/S0129054198000301
  * [37] Schweitzer, P., Schweitzer, P.: Connecting face hitting sets in planar graphs. Inf. Process. Lett. 111(1), 11–15 (2010). https://doi.org/10.1016/j.ipl.2010.10.008
  * [38] Shamir, R., Sharan, R., Tsur, D.: Cluster graph modification problems. Discret. Appl. Math. 144(1-2), 173–182 (2004). https://doi.org/10.1016/j.dam.2004.01.007
  * [39] Trotter, W.T., Harary, F.: On double and multiple interval graphs. J. Graph Theory 3(3), 205–211 (1979). https://doi.org/10.1002/jgt.3190030302
  * [40] Ueno, S., Kajitani, Y., Gotoh, S.: On the nonseparating independent set problem and feedback set problem for graphs with no vertex degree exceeding three. Discret. Math. 72(1-3), 355–360 (1988). https://doi.org/10.1016/0012-365X(88)90226-9
  * [41] Wu, H.Y., Nöllenburg, M., Viola, I.: Multi-level area balancing of clustered graphs. IEEE Trans. Vis. Comput. Graph. 28(7), 2682–2696 (2022). https://doi.org/10.1109/TVCG.2020.3038154
  * [42] Yannakakis, M.: Node-and edge-deletion NP-complete problems. In: Proc. 10th Annual ACM Symposium on Theory of Computing (STOC). pp. 253–264. ACM (1978). https://doi.org/10.1145/800133.804355