Conjecturing-Based Discovery of Patterns in Data

Conjecturing-Based Discovery of Patterns in Data111Data Ethics and
Reproducibility Note: Code and data to reproduce results are available here:
https://github.com/jpbrooks/conjecturing. COVID-19 synthetic patient data was
obtained as part of the Veterans Health Administration (VHA) Innovation
Ecosystem and precisionFDA COVID-19 Risk Factor Modeling Challenge and is used
here with permission from the Food and Drug Administration (FDA).

J. Paul Brooks Department of Information Systems, Virginia Commonwealth
University, Richmond, VA 23284<EMAIL_ADDRESS>David J. Edwards Department
of Statistical Sciences and Operations Research, Virginia Commonwealth
University, Richmond, VA 23284<EMAIL_ADDRESS>Craig E. Larson Department
of Mathematics and Applied Mathematics, Virginia Commonwealth University,
Richmond, VA 23284<EMAIL_ADDRESS>Nico Van Cleemput Department of Applied
Mathematics, Computer Science and Statistics, Ghent University, Ghent,
Belgium<EMAIL_ADDRESS>
We propose the use of a conjecturing machine that generates feature
relationships in the form of bounds involving nonlinear terms for numerical
features and boolean expressions for categorical features. The proposed
Conjecturing framework recovers known nonlinear and boolean relationships
among features from data. In both settings, true underlying relationships are
revealed. We then compare the method to a previously-proposed framework for
symbolic regression and demonstrate that it can also be used to recover
equations that are satisfied among features in a dataset. The framework is
then applied to patient-level data regarding COVID-19 outcomes to suggest
possible risk factors that are confirmed in medical literature.

automated conjecturing, computational scientific discovery, interpretable
artificial intelligence, nonlinear pattern discovery, boolean pattern
discovery

##  1 Introduction

Modern machine learning methods allow one to leverage complex relationships
present in data to generate accurate predictions but do not reveal them to the
investigator (Section 3). We propose an automated conjecturing framework for
discovering nonlinear and boolean relationships among the features in a given
dataset. Our primary goal is discovery - to provide the investigator with a
manageable number of suggested relationships to inspire future investigation
for validation.

The nonlinear relationships are produced in the form of bounds. Bounds are
useful for scientific discovery from numeric data because they 1) suggest
direct and indirect relationships among features, 2) suggest a functional form
for the relationships, and 3) can subsequently be used as boolean features
(e.g., is this bound satisfied by an observation?) for discovering more
complex boolean relationships. Whereas previous related approaches seek to
find equations for numeric data, our Conjecturing method produces bounds for
numeric data, boolean expressions for discrete data, and bounds and boolean
expressions for mixed data.

Udrescu and Tegmark (2020) proposed a system called AI Feynman that combines
deep learning with methods for symbolic regression to recover nonlinear
relationships in data. Impressively, they recover over 100 equations of
varying complexity from data. In contrast to AI Feynman, our Conjecturing
framework uses Fajtlowicz’s Dalmatian heuristic (Fajtlowicz 1995) to discover
bounds rather than equations. Further, our framework can also be applied to
categorical data to discover boolean relationships among features and already-
discovered bounds. This work represents the first application of the Dalmatian
heuristic to learning both nonlinear and boolean relationships from data. The
bounds and conditions produce interpretable yet complex relationships.

##  2 Background and Previous Related Work

In this section, we provide background on our Conjecturing framework including
examples of uses of bounds and sufficient conditions, a description of the
core algorithm, and a survey of previous related work.

###  2.1 Motivating Examples for Conjectured Bounds and Sufficient Conditions

The algorithm we use to conjecture feature relationships is an adaptation of
an algorithm that was originally designed to conjecture relationships for
mathematical objects. To illustrate the potential value of bounds and
sufficient conditions, we describe two problems and relevant results from
graph theory. This paper extends these ideas regarding bounds and sufficient
conditions to learning from data.

A graph is a collection of nodes $V$ and edges $E$ that are ordered pairs of
nodes. Consider the problem of finding bounds for the independence number of a
graph222The independence number is the largest number of nodes in a graph no
two of which are contained in an edge. The definition of independence number
is not important for this example, but only the fact that with every graph is
associated a number called the “independence number”.. It is well known that
the linear programming (LP) relaxation of an appropriate integer program
provides a upper bound on the independence number (Schrijver 2003). The Lovász
$\vartheta$ number of a graph also provides an upper bound that is known to be
no larger than the LP relaxation bound for any graph (Lovász 1979). Therefore,
the Lovász $\vartheta$ bound dominates the LP relaxation bound, and such
relationships are commonly pursued. However, relationships among bounds can be
more nuanced. Consider a third bound on the independence number due to Haemers
(1979). For some graphs it is a stronger bound than Lovász $\vartheta$ while
on other graphs it is a weaker bound; for some graphs, Lovász $\vartheta$ is a
sharp bound and Haemers’s bound is not while for other graphs, Haemers’s bound
is a sharp bound and Lovász $\vartheta$ is not. It remains an open question
whether there are a “small” number of bounds where the largest for value for
any graph would provide a sharp bound on the independence number. In this
paper, we describe a computational approach to discover bounds among numeric
features in a dataset. As with the independence number, collections of bounds
can provide valuable insight into relationships for the system from which the
data was collected.

Now consider the problem of finding sufficient conditions for a graph to be
Hamiltonian333A Hamiltonian graph is a graph with a spanning cycle (West
2001). The definition of Hamiltonian is not important for this example, but
only the fact that any graph either is or is not Hamiltonian.. Chvátal (1972)
proved that for a graph $G$ with certain conditions on the vertex degrees, $G$
is Hamiltonian. Also, Chvátal and Erdös (1972) proved that if a graph
satisfies a connectivity condition, then it is Hamiltonian. These are two
conditions that are sufficient for a graph to be Hamiltonian, but neither
implies the other. Some graphs satisfy both conditions, some graphs satisfy
one condition, and some graphs satisfy neither condition. The existence and
discovery of a (small) set of sufficient conditions that characterize all
Hamiltonian graphs remains an open area of research. The pursuit of sufficient
conditions of graph properties such as Hamiltonicity mirrors that of bounds
(Larson and Van Cleemput 2017). In the context of learning from data, we show
how categorical data, together with bounds discovered among numeric features,
can be used as input to a computational approach for generating sufficient
conditions for a property of interest.

###  2.2 The Dalmatian Heuristic

Our Conjecturing framework is based on an implementation of Fajtlowicz’s
Dalmatian heuristic (Fajtlowicz 1995, Larson and Van Cleemput 2017). The
heuristic was originally implemented in Graffiti (Fajtlowicz 1995) which was
the first program to produce research conjectures that led to new mathematical
theory. The program produces statements that are relations between
mathematical invariants which are numerical attributes of examples. Recent
implementations of the Dalmatian heuristic have been applied to the discovery
of relationships for graphs (Larson and Van Cleemput 2016) and game strategies
(Bradford et al. 2020). The heuristic was adapted to work with properties
which are boolean attributes of examples by Larson and Van Cleemput (2017). We
built our framework using a more recent implementation of the Dalmatian
heuristic available here: http://nvcleemp.github.io/conjecturing/.

We now describe invariant conjecturing using Fajtolwicz’s Dalmatian heuristic.
The inputs include the following. Let $E$ be a set of examples of a given type
(e.g., graphs or data observations). Let
$A=\\{\alpha_{1},\alpha_{2},\ldots,\alpha_{m}\\}$ be real number invariants.
In this work the examples are $n$ data observations and the invariants are $m$
numeric features. The real-numbered value of example $i$ for invariant
$\alpha_{j}$ is $\alpha_{j}(i)=x_{ij}$ for $i=1,\ldots,n$ and $j=1,\ldots,m$.
Let $O$ be a collection of unary operators and binary operators. Examples of
unary operators include adding 1, squaring, square-rooting, and division by 2.
Binary operators include addition, multiplication, and subtraction. Let
$\alpha^{*}\in A$ be the invariant for which upper and lower bounds are of
interest, and let $\alpha^{*}(i)$ be the value of the invariant of interest
for example $i$.

The aim is to generate conjectured bounds that are true for any realization of
input examples $E$. The Dalmatian heuristic provides criteria for generating
conjectured bounds that are the best for $E$. Algorithm 1 provides a way to
generate expressions of increasing complexity, apply the heuristic, and store
conjectures. The complexity of an expression is the number of nodes in the
corresponding expression tree (Figure 1) and is the sum of the number of
invariants, number of unary operators, and number of binary operators. The
algorithm proceeds by generating unlabeled trees and then labeling the nodes
with operators and invariants. Expressions satisfying the Dalmatian heuristic
conditions are retained as conjectures $\mathcal{C}$. For a conjecture $c$,
let $c(i)$ be the conjectured bound for example $i$.

With examples $E$, invariants $A$, operators $O$, invariant of interest
$\alpha^{*}$, an upper limit on the proportion of missing values allowed for
an invariant $skips$, and a direction indicating if the algorithm will produce
upper or lower bounds ($UPPER$ or $LOWER$), procedure Conjecturing-INV is
called (Algorithm 1, line 1). The number of unary nodes $u$ and binary nodes
$b$ of an expression tree are initialized to zero and the conjectures
$\mathcal{C}$ is initialized to the empty set (Algorithm 1, line 2). Line 3 of
Algorithm 1 refers to the stopping criteria of the expression generator. For
invariant conjecturing for upper bounds, if the minimum conjectured bound is
tight for each example (i.e., $\min_{c\in\mathcal{C}}c(i)=\alpha^{*}(i)$ for
$i\in E$), then the expression generator is stopped. Otherwise, expression
generation continues until a time limit is reached. Line 4 calls a procedure
to generate a tree, the branching nodes of which will be operators and the
leaf nodes are invariants. Lines 5-11 help to enumerate every tree where each
vertex connected to a leaf node has degree one or two. These branching nodes
will correspond to unary or binary operators, respectively, when the tree is
labeled. The leaf nodes will correspond to invariants. Unlabeled trees are
grown recursively and then the nodes are labeled with operators and
invariants.

The procedure generateTree (Algorithm 1, line 4) creates a new tree with a
single node, then calls procedure generateTreeRec to add new nodes until there
are $u$ unary nodes and $b$ binary nodes. The procedure generateTreeRec
(Algorithm 1, line 19) either calls generateLabeledTree to apply labels by
assigning invariants to leaf nodes and operators to branching nodes to
generate an expression (Algorithm 1, line 21), or adds nodes to grow the tree
(lines 23-32).

The procedure generateLabeledTree (Algorithm 1, line 35) takes as input a tree
with $u$ unary nodes and $b$ binary nodes. Line 36 orders the children nodes
in the order left to right appear before their parent. Then line 37 creates a
set of labeled trees. The leaf nodes are labeled with invariants and the
branching nodes are labeled with operators. Invariants with more than $skips$
missing values among examples are not used for labeling. For the commutative
binary operators, the left child is larger than the right if the left has more
nodes. If the number of nodes is equal, we use the lexicographically largest
string of labels. Since the suffix order guarantees that all subtrees are
fully labeled before their parent is labeled, this is an unambiguous
definition. Examples of labeled expression trees are given in Figure 1.

$+$$\div$$x_{1}$$x_{2}$$x_{3}$
---
(a)
$\div$$\times$pow$\times$$k$$m_{1}$$m_{2}$$r$$2$
(b)
Figure 1: Expression trees for (a) an upper bound on square footage
$x_{1}/x_{2}+x_{3}$ where $x_{1}$ is 300K, $x_{2}$ is pricePerSquareFoot, and
$x_{3}$ is bathrooms and (b) gravitational force $km_{1}m_{2}/r^{2}$.

Lines 39 and 41 are the Dalmatian heuristic. A conjectured upper bound $c$ is
only retained in the database of conjectures $\mathcal{C}$ if the bound passes
the following two tests:

  1. 1.

(Truth test). The candidate conjecture $\alpha^{*}(i)\leq c(i)$ is true for
all examples $i\in E$, and

  2. 2.

(Non-dominance test.) There is an example $i$ where
$c(i)<\min\\{c^{\prime}(i):c^{\prime}\in\mathcal{C}\setminus\\{c\\}\\}$. That
is, the candidate conjecture would give a better bound for $\alpha^{*}(i)$
than any previously conjectured (upper) bound.

Line 41 ensures that the number of conjectures is no larger than the number of
examples; i.e., $|\mathcal{C}|\leq|E|$.

The procedure is the same for generating lower bounds with the only difference
being how the Dalmatian heuristic criteria are evaluated in lines 39 and 41.

The computational requirements of Algorithm 1 increase exponentially with the
number of invariants and with the number of operators. The computation time
increases with the number of examples because of the check in Step 39. To
facilitate generation of more candidate expressions in less time, one can use
fewer examples as input to the algorithm. To achieve additional efficiency, we
implement the following design choice. In our implementation, when a tree is
labeled, operators can be reused, but invariants cannot. We make this design
choice so that more expressions can be generated in a smaller amount of time.
In Section 5.4, we will demonstrate how this limitation can be overcome in
situations where repeating invariants is warranted.

Algorithm 1 Invariant Conjecturing.

Input: Examples $E$, Invariants $A$, operators $O$, invariant of interest
$\alpha^{*}$, invariant missing value limit $skips$, direction ($UPPER$ or
$LOWER$).
Output: Conjectured $\mathcal{C}$ in the form of conjectured bounds on the
invariant of interest $\alpha^{*}$.

1:procedure Conjecturing-INV

2: Set $u=0$, $b=0$, $\mathcal{C}=\emptyset$.

3: while not stopped do

4: generateTree($u$, $b$).

5: if $b=0$ then

6: Set $b=\lceil u/2\rceil$.

7: Set $u=0$.

8: else

9: $b--$.

10: $u+=2$.

11: end if

12: end while

13: return $\mathcal{C}$.

14:end procedure

15:procedure generateTree($u$, $b$)

16: Set tree = new tree with single node.

17: generateTreeRec(tree, $u$, $b$).

18:end procedure

19:procedure generateTreeRec(tree, $u$, $b$)

20: if number of unary nodes == $u$ and number of binary nodes == $b$ then

21: generateLabeledTree(tree).

22: else

23: for all nodes $v$ on the second-deepest level that have at most 1 child
and have no nodes at the same level to their right with at least 1 child do

24: Add child to $v$.

25: generateTreeRec(tree, $u$, $b$).

26: Remove that child from $v$.

27: end for

28: for all nodes $v$ on the deepest level do

29: Add child to $v$.

30: generateTreeRec(tree, $u$, $b$).

31: Remove that child from $v$.

32: end for

33: end if

34:end procedure

35:procedure generateLabeledTree(tree)

36: Order the nodes in a suffix order.

37: Recursively label each node in this ordered array with either an
invariant, a unary operator, or a binary operator depending on its degree. For
commutative binary operators we only label a vertex if its left child is
larger than its right child.

38: for each fully labeled tree do

39: if the corresponding bound $c$ is valid for all examples in $E$ and is not
dominated by existing bounds in $\mathcal{C}$. then

40: Set $\mathcal{C}=\mathcal{C}\cup{c}$.

41: Remove dominated conjectures from $\mathcal{C}$.

42: end if

43: end for

44:end procedure

---
(a)
(b)
Figure 2: (a) Upper bounds and (b) lower bounds generated for gravitational
force using Conjecturing-INV, the invariant version of the conjecturing
algorithm. Instances from the training data are on the $x$-axis. The gold
curve is the true value for the instances. The blue curve in (a) is the true
value without the constant of proportionality and is one of the upper bounds.
The red curve is the (a) maximum and (b) minimum of the discovered bounds.

Figure 2 displays (a) upper bounds and (b) lower bounds derived for test
instances for data generated based on a formula for gravity. The gray curves
correspond to bounds, and each must be the best on at least one training
example instance in order to be retained. More details on this experiment are
provided in Section 3.1 and Section 4.1.

Algorithm 1 can be adapted for property conjecturing with few modifications.
We now detail the differences. Let $E$ be a set of examples and let
$\Pi=\\{\pi_{1},\pi_{2},\ldots,\pi_{m}\\}$ be properties. The examples are $n$
data observations and the properties are $m$ boolean features. The truth value
of example $i$ for property $\pi_{j}$ is $\pi_{j}(i)$. Let $O$ be the
following collection of operators: NOT ($\neg$), AND ($\&$), OR ($|$), XOR
(exclusive or) ($\oplus$), and IMPLIES ($\rightarrow$). NOT is a unary
operator and the remaining operators are binary operators. Let $\pi^{*}\in\Pi$
be the property for which sufficient and/or necessary conditions are of
interest, and let $\pi^{*}(i)$ be the truth value of the property of interest
for example $i$.

The aim is to generate conjectured sufficient or necessary conditions for the
property of interest that are valid for any realization of input examples $E$.
The algorithm for property conjecturing procedure Conjecturing-PROP generates
unlabeled trees as in Algorithm 1 but then labels the nodes with operators and
properties. Logical expressions satisfying the Dalmatian heuristic conditions
are retained as conjectures $\mathcal{C}$. For a conjecture $c$, let $c(i)$ be
the conjectured truth value for example $i$.

The inputs to property conjecturing are examples $E$, properties $\Pi$,
operators $O$, a property of interest $\pi^{*}$, and a direction
($SUFFICIENT$, $NECESSARY$) indicating if the algorithm will produce
sufficient or necessary conditions for the property of interest.

The stopping criterion for property conjecturing for the case that direction
is $SUFFICIENT$ is obtaining a set of conjectures where for every example with
$\pi^{*}(i)=$ true, each example evaluates to true for at least one
conjecture. Otherwise, expression generation continues until a time limit is
reached.

The Dalmatian heuristic for property conjectures is applied as follows. A
conjectured sufficient condition $c$ is only retained in the database of
conjectures $\mathcal{C}$ if the expression passes the following two tests:

  1. 1.

(Truth test). For all examples $i\in E$ for which $c(i)$ is true, then
$\pi^{*}(i)$ is also true, and

  2. 2.

(Non-dominance test.) The number of examples $i\in E$ for which $c(i)$ is true
is not a subset of examples that evaluate to true for any previously
conjectured sufficient condition.

To generate necessary conditions for $\pi^{*}$, one can generate sufficient
conditions for $\neg\pi^{*}$ (NOT $\pi^{*}$).

---
(a)
(b)
Figure 3: Schematic of property conjecturing. In (a), conditions 1 and 2
evaluate to true for a subset of examples with the property of interest and
for no samples that do not have the property of interest. Condition 3
evaluates to true for examples with and without the property of interest, and
so it is discarded. In (b), the union of sufficient conditions covers all
examples with the property of interest.

Figure 3(a) depicts candidate conditions for examples with the property of
interest (green) and those without (red). Conditions 1 and 2 are sufficient
conditions for subsets of examples with the property of interest. Condition 3
evaluates to true for examples with and without the property of interest, and
would therefore not be retained. The goal of property conjecturing is to find
a set of sufficient conditions that evaluate to true for all examples with the
property of interest and for none of the examples without the property of
interest as illustrated in Figure 3(b).

###  2.3 Other Related Work

In this section, we explain how our work is related to previous work in
automated scientific discovery, machine learning interpretability, automated
feature engineering, and empirical model building.

Symbolic regression has been used as a tool for automated scientific
discovery. Symbolic regression is the use of genetic programming to
approximate a target function on training data and generalize to produce
predictions on new data (Nicolau and Agapitos 2021) and until the work of
Schmidt and Lipson (2009), the focus was on improving prediction accuracy by
approximating an underlying function rather than a focus on discovering true
functional relationships among features. Schmidt and Lipson (2009) extend
previous work to develop a system for discovering laws for dynamical systems
by considering relationships among derivatives. Their work led to the
development of a software Eureqa. More recently, Udrescu and Tegmark (2020)
combined a variety of strategies including dimensional analysis, symmetry
identification, neural network training, and brute-force enumeration into a
framework called AI Feynman to recover true physical functional forms from
data. Petersen et al. (2021) propose a method for deep symbolic regression
that combines reinforcement learning with a recurrent neural network model.
They compare their approach with methods based on priority queue training
proposed by (Abolafia et al. 2018) and the traditional genetic programming
approach. Our experiments include comparisons with each of these methods in
the ability to recover equations from data.

Other frameworks have been proposed to computationally generate conjectures
from data and discover scientific laws. Data smashing is introduced by
Chattopadhyay and Lipson (2014) as a method for computing dissimilarities from
streams of data (e.g., electroencephalogram data) to aid in revealing
relationships among observations. Jantzen (2020) proposes an algorithm with
the similar purpose of detecting types of dynamical systems called dynamical
kinds. Subsequently, these kinds “are then targets for law-like
generalization” (Jantzen 2020). While Jantzen’s work provides a method for
discovering the kinds, it does not suggest how to recover the “laws”. It is
these relationships that we aim to discover with the Conjecturing framework.

Our work is distinguished from these previous works in that 1) we focus on
generating bounds for invariants that serve as hypotheses for the investigator
rather than recover true functional forms or generate accurate predictions, 2)
our invariant conjecturing algorithm is paired with a property conjecturing
algorithm for discovering both nonlinear bounds and boolean relationships, 3)
our framework is designed for a given static observational dataset rather than
on discovering laws for dynamical systems, and 4) rather than a stochastic
search over the space of functional forms, our Conjecturing system leverages
sophisticated techniques for enumerating expressions of increasing complexity
(described in Larson and Van Cleemput (2016) for “noiseless” data involving
mathematical objects such as graphs). In our system, the human remains “in the
loop” to evaluate the plausibility of suggested bounds and conditions.

Brunton et al. (2016) introduce SINDy for combining sparse regression and
expert knowledge to developing models of dynamical systems. We adopt a similar
approach to incorporating prior knowledge in Section 5.4 for the Nguyen
benchmark suite (Nguyen et al. 2011) where we provide the Conjecturing
framework with candidate non-linear functions as building blocks. Unlike our
Conjecturing framework, theirs is designed for recovering equations governing
dynamical systems rather than bounds and theirs is not capable of recovering
boolean relationships.

Langley (2019) provides a review of past efforts in computational scientific
discovery. Several frameworks have their origins in analyzing mass
spectroscopy and other electrochemical data. Bacon (Langely et al. 1987) is a
general framework for scientific discovery based on suggesting and executing a
series of designed experiments. Tallorin et al. (2018) proposed a method
called POOL that uses Bayesian optimization and machine learning in an
iterative fashion for experiments to discover peptide substrates for enzymes.
Bacon and POOL both make recommendations regarding additional data to collect
while our system assumes that a fixed dataset is provided that may or may not
be the result of a designed experiment.

Precise definitions of “explainability” and “interpretability” are still being
developed (Vilone and Longo 2020, Lu et al. 2019, Fürnkranz et al. 2020) as
research in the area has rapidly accelerated. According to the convention of
Rudin (2019), explainability is concerned with post-hoc analyses of black box
models to create simple explanations of model behavior. Motivated by observed
accuracies of deep learning models, work in this area includes identifying
important features for prediction, building simple local models, conducting
sensitivity analyses, and deriving prototype examples (Samek and Müller 2019,
Elton 2020). Tsang et al. (2018a, b, 2020) develop neural network frameworks
for identifying sets of features for which there is an interaction \- a non-
additive relationship among predictive features that influence a response
value. These methods provide explainability in that they identify sets of
features that interact, but the framework is not designed to reveal the
functional form of the nonlinear interaction.

Rudin (2019) advocates the development of interpretable models where the
mechanisms for predictions are simple relationships that are readily apparent
to the investigator. Much of the recent work in this area is in the
development of decision rules (e.g., (Hammer and Bonates 2006, Dash et al.
2018, Bellomarini et al. 2020)) or decision lists and trees (e.g., (Wang and
Rudin 2015, Wang et al. 2017, Rudin and Ertekin 2018, Bertsimas and Dunn 2017,
Verwer and Zhang 2019, Blanquero et al. 2021, Aghaei et al. 2021, Akyüz and
Birbil 2021)). Different from these works, our Conjecturing framework
automates the discovery of nonlinear features. In addition, as with work on
decision rules in general, our framework can combine the discrete features in
data with the discovered nonlinear features to discover a potentially richer
set of boolean relationships when compared to optimization-based trees and
decision lists.

Khurana et al. (2018) propose a system that leverages reinforcement learning
to search expression trees for predictive features. ExploreKit (Katz et al.
2016) is a framework for automatic feature engineering that combines features
using basic arithmetic operations and then uses machine learning to rank and
select those with high predictive ability. The Data Science Machine (Kanter
and Veeramachaneni 2015) automatically generates features for entities in
relational databases with possible dependencies between tables followed by
singular value decomposition. In none of these works is model transparency
evaluated but rather only model performance. An important distinction of our
work from these is that they focused on improving prediction accuracy,
sometimes at the expense of understandable features, and not on scientific
discovery.

Traditional statistical methods for empirical model building (e.g. regression
analysis) tend to focus on first- and second-order polynomial models;
interaction terms up to a certain degree are often included. Empirical models
are intended to provide adequate prediction performance while also providing a
simple assessment of feature importance via model coefficients. Techniques
such as all-subsets, stepwise selection, and regularization methods (e.g.,
LASSO (Tibshirani 1996)) are commonly used to perform feature selection over
model spaces of increasing complexity. However, domain knowledge is typically
required for reciprocal or non-polynomial relationships. Our Conjecturing
framework provides a search over a much broader class of nonlinear functions.

##  3 Two Motivating Examples

In this section, we describe two datasets where a “typical” knowledge
discovery workflow fails to reveal important relationships among features.

Research on machine learning does, of course, lead to conjectured
relationships between variables which are in turn used to make predictions of
one or more variables in terms of others. A trained neural net, for instance,
can be viewed as a black box representing a function which produces an output
for every input in its domain. These functions are complex and of a different
character than classical scientific laws: in particular, there is little hope
of deriving these functions or relationships from simpler existing laws. Our
Conjecturing framework aims to help fill this gap in current capabilities.

###  3.1 Discovering Gravity

In this example, a numeric invariant of interest is determined by a more
complex nonlinear relationship with three numeric predictors. Consider
measurements including the masses of two objects $m_{1}$ and $m_{2}$, their
distance $r$, and the gravitational force between them $F$. The goal is to
recover the dependence of $F$ on $m_{1}$, $m_{2}$, and $r$, or

$F=k\frac{m_{1}m_{2}}{r^{2}},$

where $k$ is the gravitational constant. Following the demonstration by
Langely et al. (1987), we create a fictional dataset using a predefined value
for $k$ that is a random number between 0 and 1. For our illustrative example,
we generated 1,000 training data points and 1,000 test data points with
$k=0.057098$. Values for $m_{1}$, $m_{2}$, and $r$ are samples from
Uniform(1,100000) distributions, and $F$ is calculated for each sample with no
noise.

A linear regression model will fail to capture the nonlinear interaction of
the variables. Off-the-shelf machine learning methods such as random forests
and neural networks can leverage the nonlinear relationship in the data but
cannot present the relationship to the investigator. In the next section, we
propose a framework for producing bounds on $F$ that are functions of the
other features.

###  3.2 Discovering an Interaction in Real Estate Valuation Data

The second example is a case where a boolean variable of interest is
completely determined by the product of two numeric features in the dataset;
i.e., the second-order interaction term completely defines the relationship.

Consider a dataset on residential real estate properties for sale obtained
from https://www.redfin.com. The goal is to predict whether a home with given
feature values has a list price above or below $300,000.

This dataset includes both the price per square foot and total square footage
along with eight additional features such as the number of bathrooms and
bedrooms. The property of interest (above vs. below) can be determined by
multiplying the price per square foot by square footage and setting a
threshold. Thus, the interaction of price per square foot and square footage,
hereafter called the active interaction, completely describes the relationship
between the predictors and response. Data are partitioned into a training
dataset with 1,000 houses and a test dataset with 30,156 houses. In the next
section, we leverage our framework for invariant bounds and then extend it to
produce boolean relationships to recover the active interaction term and how
it determines class membership.

##  4 A Conjecturing Framework for Discovering Patterns in Data

We now describe a framework that leverages a conjecturing algorithm to
discover nonlinear and boolean feature relationships in data. All experiments
were run on a computer with an Intel i7-2600 CPU @ 3.4GHz and 16 GB RAM.

###  4.1 Conjecturing for Nonlinear Relationships

The invariant version of the conjecturing method (procedure Conjecturing-INV)
can be used for discovering nonlinear relationships in data. Invariant
conjectures are generated that provide upper and lower bounds on the invariant
of interest. These conjectures are the nonlinear functions that can be used as
new features and/or as a complete model for the system.

For the gravity case from Section 3.1, the invariants are
$A=\\{F,m_{1},m_{2},r\\}$ and the invariant of interest is the force $F$. The
examples $E$ are the observations in the data.

The Conjecturing framework is not designed to recover constants such as the
gravitational constant $k$. In general, for a functional relationship with a
constant $k$ such that $0<k<1$, the expression without the constant provides a
lower bound for the response and the reciprocal expression provides an upper
bound. In cases where the constant is not between 0 and 1, the converse is
true. For this example, $F=km_{1}m_{2}/r^{2}$ with $0<k<1$ and so
$r^{2}/m_{1}m_{2}\leq F\leq m_{1}m_{2}/r^{2}$. The Conjecturing framework can
potentially recover the bounds $F\leq m_{1}m_{2}/r^{2}$ or $F\geq
r^{2}/m_{1}m_{2}$.

For our example, Conjecturing-INV returns 19 upper bounds and 24 lower bounds
for $F$. Among the upper bounds is

$F\leq m_{1}m_{2}/r^{2},$

which approximates the true gravity relationship used to generate the data.
The bound does not include the constant $k$. The lower bound of
$r^{2}/m_{1}m_{2}$ is not recovered. Other bounds generated by Conjecturing-
INV include

$\displaystyle F$ $\displaystyle\leq 2m_{2}/\sqrt{r},$ (1) $\displaystyle F$
$\displaystyle\leq 2|m_{1}-m_{2}|,$ (2) $\displaystyle F$ $\displaystyle\geq
8m_{2}/r^{2},$ (3) $\displaystyle F$ $\displaystyle\geq-1/(r-2m_{2}).$ (4)

Eight of the upper bounds and 15 of lower bounds for $F$ are depicted in
Figure 2. The upper bound $m_{1}m_{2}/r^{2}$ in Figure 2(a) is blue, while the
true value $km_{1}m_{2}/r^{2}$ is gold.

As the primary goal of our approach is discovery, the bounds produced are
suggestions that require further validation. We consider it a success that the
relationship $F\propto m_{1}m_{2}/r^{2}$ is included in one of the bounds.
Among the other bounds produced, we see that true relationships $F\propto
m_{1}$, $F\propto m_{2}$, and $F\propto 1/r^{2}$ are all suggested, along with
false relationships $F\propto 1/\sqrt{r}$ and $F\propto 1/r$. Follow up
investigations can be used to inspect these relationships and potentially
recover the gravitational constant. The lower bound that is not recovered,
while valid, does not reflect the true proportionality relationships which may
be why it is dominated by other bounds. An approximation of the gravitational
constant of 0.057 could be represented as $(+1+1+1+1+1)\times
10^{-1-1}+(+1+1+1+1+1)\times 10^{-1-1-1}$ which by itself has complexity 19.
The expression $m_{1}m_{2}/r^{2}$ has complexity 6.

In this example, the Conjecturing framework recovers the true nonlinear
relationship up to a constant of proportionality along with 44 additional
suggested bounds. Therefore, isolating a single true bound, in the case where
the bound is unknown, can require additional analysis and/or experiments. The
additional bounds can provide potential insight into feature interactions.

###  4.2 Conjecturing for Nonlinear and Boolean Relationships with Mixed Data

Our Conjecturing framework for mixed data leverages the invariant version
(procedure Conjecturing-INV) and the property version (procedure Conjecturing-
PROP) of the conjecturing algorithm. For mixed data, we propose a framework to
produce conjectures of nonlinear and boolean patterns. These conjectures can
capture complex patterns while maintaining interpretability.

We assume that we are given a dataset with numeric features $N$, boolean
features $B$, and a categorical feature of interest with levels $\mathcal{Y}$.
Note that a categorical feature with more than two levels can be converted to
a series of boolean features. Let $\pi_{y}$ be the property that an
observation has value $y$, for $y\in\mathcal{Y}$.

For each level $y\in\mathcal{Y}$, the algorithm discovers bounds for the
numeric features that are satisfied by each observation in the class
(Algorithm 2, Lines 4-14). These inequalities are converted to properties of
the form “if the inequality is satisfied, then true; false, otherwise”
(Algorithm 2, Line 12). These new properties are combined with the original
boolean features in the data (Algorithm 2, Line 13). The properties from
across all classes are pooled together and the observations belonging to all
classes are pooled together as examples and then, for each level
$y\in\mathcal{Y}$, the property version of conjecturing is applied to discover
sufficient conditions for $\pi_{y}$ (Algorithm 2, Lines 15-20).

We now provide further details on Algorithm 2 using the real estate valuation
case from Section 3.2 as an illustrative example. First, we convert the
categorical feature propertyType into boolean features condo, mobileHome,
singleFamily, townhouse, multiFamily2-4Unit, multifFamily5PlusUnit, and Other.
We also add a feature that is a constant value of 300,000 for each observation
because it is the price cutoff and call it 300K. The resulting $18$ features
are partitioned into numeric features $N=\\{$bedrooms, bathrooms,
squareFootage, lotSize, yearBuilt, daysOnMarket, pricePerSquareFoot,
hoaPerMonth, latitude, longitude, 300K$\\}$ (Algorithm 2, Line 2) and boolean
features $B=\\{$condo, mobileHome, singleFamily, townhouse,
multiFamily2-4Unit, multifFamily5PlusUnit, Other$\\}$ (Algorithm 2, Line 3).

In our training set, there are 1,000 observations that are used as examples.
For each value of the property of interest, $\\{\textrm{\it below},\textrm{\it
above}\\}$, the corresponding observations serve as the examples (Algorithm 2,
Line 6). For each numeric feature, upper and lower bounds on that feature are
found that are functions of the other numeric features (Algorithm 2, Lines
8-9). These are found by applying the invariant relations version of the
conjecturing method (Conjecturing-INV). For houses with property below, there
are 1,280 bounds derived. Included are plausible relations concerning house
features that are seemingly irrelevant to the classification task such as

bathrooms $\displaystyle\leq 2\times\textit{bedrooms}$ (5) bedrooms
$\displaystyle\geq\textit{bathrooms}-1$ (6) lotSize
$\displaystyle\geq(\textit{squareFootage}-\textit{yearBuilt})\times\textit{bedrooms}$
(7)

Also included are less-interpretable bounds such as:

yearBuilt
$\displaystyle\geq\textit{hoaPerMonth}\times\log(10)/\log(2\times\textit{daysOnMarket})$
(8) daysOnMarket $\displaystyle\leq e^{e^{\sqrt{2\times\textit{lotSize}}}}$
(9) hoaPerMonth $\displaystyle\leq
10^{2\times\textit{bathrooms}}+\textit{squareFootage}$ (10)

There were also several bounds discovered that are close approximations of the
relationship present in the active interaction term, including

squareFootage
$\displaystyle\leq\textit{300K}/\textit{pricePerSquareFoot}+\textit{bathrooms}$
(11) squareFootage
$\displaystyle\leq\textit{300K}/\textit{pricePerSquareFoot}+\textit{bedrooms}$
(12) squareFootage
$\displaystyle\leq\textit{300K}/\textit{pricePerSquareFoot}+\textit{daysOnMarket}$
(13) squareFootage
$\displaystyle\leq\textit{300K}/(\textit{pricePerSquareFoot}-1)-1$ (14)
pricePerSquareFoot
$\displaystyle\leq-\textit{300K}/(\textit{bedrooms}-\textit{squareFootage})$
(15) pricePerSquareFoot
$\displaystyle\leq\left\lceil\textit{300K}/\textit{squareFootage}\right\rceil$
(16) 300K
$\displaystyle\geq-(\textit{bathrooms}-\textit{squareFootage})\times\textit{pricePerSquareFoot}$
(17)

For houses with property above, there are 1,457 bounds derived including a mix
of simple relations and less intuitive relations. Also included are the
following three relations that are nearly identical to the active interaction
relation:

squareFootage $\displaystyle\geq\textit{300K}/(\textit{pricePerSquareFoot}+1)$
(18) pricePerSquareFoot
$\displaystyle\geq\textit{300K}/\textit{squareFootage}+1$ (19) 300K
$\displaystyle\leq(\textit{pricePerSquareFoot}+1)\times\textit{squareFootage}$
(20)

The resulting invariant relations are pooled together (Algorithm 2, Line 10).
The invariant relations are encoded as properties (Algorithm 2, Line 12). The
original binary features from the data are also encoded as properties for a
total of $1,280+1,457+7=2,744$ properties. Examples of encoded properties from
the invariant relations are:

bathrooms
$\displaystyle\stackrel{{\scriptstyle?}}{{\leq}}2\times\textit{bedrooms}$ (21)
$\displaystyle(\textit{yearBuilt}$
$\displaystyle\stackrel{{\scriptstyle?}}{{\geq}}\textit{hoaPerMonth}\times\log(10)/\log(2\times\textit{daysOnMarket}))$
(22) $\displaystyle(\textit{squareFootage}$
$\displaystyle\stackrel{{\scriptstyle?}}{{\leq}}\textit{300K}/\textit{pricePerSquareFoot}+\textit{bathrooms})$
(23) $\displaystyle(\textit{squareFootage}$
$\displaystyle\stackrel{{\scriptstyle?}}{{\geq}}\textit{300K}/(\textit{pricePerSquareFoot}+1)).$
(24)

These properties can be used as boolean features that indicate whether a
nonlinear relationship among numeric features is satisfied for an observation.

The properties generated for each level $\\{\textrm{\it below},\textrm{\it
above}\\}$ are collected in a set $\Pi$ along with $\pi_{y}$ and the seven
original boolean features (Algorithm 2, Line 13).

For each level $\\{\textrm{\it below},\textrm{\it above}\\}$, apply the
property version of conjecturing to the properties $\Pi$ with the training
data observations serving as the examples $E$ and level as the property of
interest (Algorithm 2, Lines 15-20). The result is a set of properties that
are sufficient conditions for the levels.

Conjecturing-PROP returns only two properties. They both approximate the
underlying active interaction.

bathrooms
$\displaystyle\geq-\textit{300K}/\textit{pricePerSquareFoot}+\textit{squareFootage}$
$\displaystyle\rightarrow\textit{isBelow}$ (25) squareFootage
$\displaystyle\geq(\textit{300K}+1)/(\textit{pricePerSquareFoot}-1)$
$\displaystyle\rightarrow\textit{isAbove}$ (26)

An inspection of the data reveals that for some of the houses, there is some
rounding error when comparing the price to the square footage multiplied by
the price per square foot. The conjecturing algorithm compensates by using
invariants as error terms. In the first property, the error term is
$\textit{bathrooms}\times\textit{pricePerSquareFoot}$. In the second property,
the error term is $\textit{squareFootage}+1$.

When these properties are applied as classification rules for predicting
whether a house will be above or below $300,000, they produce no error on the
training data. The first property misclassifies 37 of 30,156 houses in the
test data for an accuracy of 0.999. The second property misclassifies 26
houses. The misclassified houses are due to rounding error and miscoding of
data. For example, one house in the test data is listed as having 31,248
bathrooms and another is listed as having a price of $459.

Despite the noise and rounding error in the data, the Conjecturing framework
was able to recover the active interaction term and these properties can be
used as features for classifiers with near-perfect accuracy.

Algorithm 2 Conjecturing framework for nonlinear and boolean relationships
with mixed data

Input: Data observations $\\{1,\ldots,n\\}$ with numeric features $N$, boolean
features $B$, and a categorical feature of interest with levels $\mathcal{Y}$;
a set of invariant operators $O$ and a set of property operators $P$.
Output: A set of conjectured properties $\mathcal{P}$.

1:Set $\mathcal{P}=\emptyset$. /* Initialize properties set. */

2:Set $A=\\{\alpha_{j}:j\in N\\}$. /* Define the set of invariants to be the
original numeric features in the data. */

3:Set $\Pi=\\{\pi_{j}:j\in B\\}$. /* Define the set of properties to be the
original boolean features in the data. */

4:for $y\in\mathcal{Y}$ do /* Loop on the levels of the categorical feature of
interest. */

5: Set $\mathcal{R}=\emptyset$. /* Initialize invariant relations set. */

6: Set $E=\\{i:\pi_{y}\\}$. /* Define the set of observations with level $y$
as the examples. */

7: for $j\in N$ do /* Loop on original numeric features. */

8: Set $R_{U}=\mbox{{Conjecutring-INV}}(E,A,O,\alpha_{j},UPPER)$ /* Submit
examples, invariants, and the invariant of interest to the invariant version
of Conjecturing for upper bounds. */

9: Set $R_{L}=\mbox{{Conjecturing-INV}}(E,A,O,\alpha_{j},LOWER)$ /* Submit
examples, invariants, and the invariant of interest to the invariant version
of Conjecturing for lower bounds. */

10: Set $\mathcal{R}=\mathcal{R}\cup R_{U}\cup R_{L}$.

11: end for

12: Convert the new invariant relations $\mathcal{R}$ into properties
$\Pi_{\mathcal{R}}$.

13: Set $\Pi=\Pi\cup\pi_{y}\cup\Pi_{\mathcal{R}}$. /* Define the set of
properties to be the original boolean features, the level $y$, and the
invariant relations properties. */

14:end for

15:for $y\in\mathcal{Y}$ do /* Loop again on the levels of the categorical
feature of interest. */

16: Set $E=\\{1,\ldots,n\\}$. /* Use all examples. */

17: Set $\mathcal{P}_{S}=\mbox{{Conjecturing-
PROP}}(E,\Pi,P,\pi_{y},SUFFICIENT)$. /* Submit examples, properties, and the
level $y$ as the property of interest to the property version of Conjecturing
for sufficient conditions. */

18: Set $\mathcal{P}_{N}=\mbox{{Conjecturing-
PROP}}(E,\Pi,P,\pi_{y},NECESSARY)$. /* Submit examples, properties, and the
level $y$ as the property of interest to the property version of Conjecturing
for necessary conditions. */

19: Set $\mathcal{P}=\mathcal{P}\cup\mathcal{P}_{S}\cup\mathcal{P}_{N}$.

20:end for

21:return $\mathcal{P}$.

Hereafter, we use “Conjecturing framework” to imply:

  1. 1.

In the case that all features are numeric, apply procedure Conjecturing-INV.

  2. 2.

In the case that all features are categorical, convert the features to a
series of properties (boolean features) and apply procedure Conjecturing-PROP.

  3. 3.

In the case of mixed data, apply Algorithm 2.

If there is no invariant of interest or no property of interest, each
invariant and/or property can serve as the invariant/property of interest in
turn, and conjectures can be generated for each.

##  5 Additional Computational Experiments

###  5.1 Sensitivity to the Number of Features

To investigate the impact of the number of features on the performance on the
Conjecturing framework, we conduct experiments adding noisy features to the
gravity example. We use the same experimental setup as described in Section
4.1 including a time limit of five seconds. For $k=0,\ldots,10$, we add $k$
noise invariants generated from a standard normal distribution and check 1)
whether Conjecturing-INV recovers $m_{1}m_{2}/r^{2}$, 2) the number of
conjectures produced, 3) the number of expressions evaluated, and 4) the
number of valid expressions produced. Valid expressions are bounds that are
valid for all training examples.

For $k=0,\ldots,6$, Conjecturing-INV recovers $m_{1}m_{2}/r^{2}$ as a
conjectured upper bound and for $k=7,\ldots,10$ the bound is not recovered
because of the additional noise invariants. The number of original invariants
is five including the force $F$, so in this experiment we can add more than
100% additional invariants and still recover the true proportionality
relationship $F\propto m_{1}m_{2}/r^{2}$.

Figure 4 contains plots of the number of conjectures produced, number of
expressions evaluated, and the number of valid expressions produced within the
five second time limit. The number of conjectures produced and expressions
evaluated increases as the number of columns increases and tends to be larger
for lower bounds than upper bounds. The number of conjectures prodcued ranges
from 22 to 174. The number of valid expressions fluctuates between 75,000 and
135,000 and there is no discernible pattern effect of the number of noise
invariants $k$. As $k$ increases, there are more invariants available and the
number of low-complexity expressions increases exponentially as does the
number of low-complexity expressions comprised of the noise invariants. Low-
complexity expressions can be generated and checked more quickly which is why
the number of expressions evaluated increases with $k$. The number of
expressions evaluated in five seconds ranges between about 300,000 and 1.4
million.

Figure 4: The number of conjectures produced, expressions tested, and valid
experessions produced as a function of the number of noise invariants added to
the gravity experiment.

###  5.2 Sensitivity to Training Examples

To investigate the effect of different subsets of training examples on the
ability of the Conjecturing framework to recover true relationships, we apply
the framework to the real estate experiment with 10 random samples of 1,000
training examples. We use the same experimental setup as described in Section
2 including a time limit of five seconds.

Table 1 contains the conjectures produced by the Conjecturing framework for
each of the 10 replications. As in the experiment in Section 4.2, the
Conjecturing framework makes use of invariants and operators to compensate for
rounding error. The invariants employed as tolerances are bathrooms and
longitude. The framework also employs operators $+1$, $-1$, and
$\lceil\cdot\rceil$ to account for deviations from the underlying active
interaction. Each of the bounds can be rewritten in terms of squareFootage
$\times$ pricePerSquareFoot $-300K$ plus or minus a small error term
containing at most one additional invariant. Therefore, we see that in this
instance, the method is not sensitive to the choice of training examples.
Further experiments are needed to understand how well the Conjecturing
framework can recover underlying relationships in the presence of different
kinds of noise. This is the subject of future work.

Table 1: Sufficient Conditions Produced by Conjecturing for Different Training Set Samples for the Real Estate Valuation Data Sample | Conjectures Produced
---|---
1 | $bathrooms\geq-300K/pricePerSquareFoot+squareFootage\rightarrow isBelow$
| $squareFootage\geq(300K+1)/(pricePerSquareFoot-1))\rightarrow isAbove$
2 | $squareFootage\leq\lceil 300K/pricePerSquareFoot\rceil\rightarrow isBelow$
| $squareFootage\geq(300K+1)/(pricePerSquareFoot-1)\rightarrow isAbove$
3 | $bathrooms\geq-300K/pricePerSquareFoot+squareFootage\rightarrow isBelow$
| $squareFootage\geq(300K+1)/(pricePerSquareFoot-1))\rightarrow isAbove$
4 | $bathrooms\geq-300K/pricePerSquareFoot+squareFootage\rightarrow isBelow$
| $squareFootage\geq(300K+1)/(pricePerSquareFoot-1))\rightarrow isAbove$
5 | $bathrooms\geq-300K/pricePerSquareFoot+squareFootage\rightarrow isBelow$
| $squareFootage\geq(300K+1)/(pricePerSquareFoot-1))\rightarrow isAbove$
6 | $squareFootage\leq\lceil 300K/pricePerSquareFoot\rceil\rightarrow isBelow$
| $squareFootage\geq(300K+1)/(pricePerSquareFoot-1)\rightarrow isAbove$
7 | $squareFootage\leq(300K+longitude)/pricePerSquareFoot\rightarrow isBelow$
| $squareFootage\geq(300K+1)/(pricePerSquareFoot-1))\rightarrow isAbove$
8 | $bathrooms\geq-300K/pricePerSquareFoot+squareFootage\rightarrow isBelow$
| $squareFootage\geq(300K+1)/(pricePerSquareFoot-1))\rightarrow isAbove$
9 | $squareFootage\leq(300K-daysOnMarket)/pricePerSquareFoot\rightarrow isBelow$
| $squareFootage\geq(300K+1)/(pricePerSquareFoot-1)\rightarrow isAbove$
10 | $bathrooms\geq-300K/pricePerSquareFoot+squareFootage\rightarrow isBelow$
| $squareFootage\geq(300K+1)/(pricePerSquareFoot-1))\rightarrow isAbove$

###  5.3 Comparison to AI Feynman (Udrescu and Tegmark 2020)

In this section, we compare the ability of Conjecturing to recover equations
from datasets used by Udrescu and Tegmark (2020) with their algorithm AI
Feynman. We then apply the implementation of AI Feynman to the gravity and
real estate datasets described in Section 3. We note that the primary goal of
Conjecturing is for discovery of nonlinear and boolean relationships while the
primary goal of AI Feynman is recovery of equations.

Performance on Feynman Equations. We apply Conjecturing to the first 10
equations listed in Table 4 of (Udrescu and Tegmark 2020) to draw comparisons
based on solution time and noise tolerance. We used the data published by the
authors here: https://space.mit.edu/home/tegmark/aifeynman.html. As in
(Udrescu and Tegmark 2020), for each instance we apply Conjecturing with three
subsets of operators in turn: $\\{+,-,\times,\div,+1,-1,^{2},\sqrt{}\\}$,
$\\{+,-,\times,\div,+1,-1,^{2},\sqrt{},\sin,\cos,\ln,^{-1},e\\}$,
$\\{+,-,\times,\div,+1,-1,^{2},\sqrt{},\sin,\cos,\ln,^{-1},e,|\cdot|,\sin^{-1},\tan^{-1}\\}$.
For instances where an equation includes the constant $\pi$, we include $\pi$
as a constant invariant. For each instance, we use the first 10 samples in
each dataset and run Conjecturing for 10, 100, and 1,000 seconds. We also run
Conjecturing for the noise tolerance of and time required by AI Feynman to
recover the equations as reported in Table 4 of (Udrescu and Tegmark 2020).

Tables 2 and 6 contain the results of applying Conjecturing to the datasets.
Conjecturing produces bounds that match the equation for four of the ten
instances. Conjecturing finds a match for all equations with complexity 10 or
less and is unable to find a match for equations with higher complexity.
Udrescu and Tegmark (2020) report that AI Feynman resolves all of the
equations while Eureqa (Schmidt and Lipson 2009) resolves four of the ten
equations, three of which are different from those found by Conjecturing.
These results indicate that Conjecturing is well suited for recovering
equations of complexity 10 or less within 1,000 seconds and sometimes within
much shorter times. Higher-complexity formulas with more invariants require
additional time.

Equations I.6.2 and I.6.2.b in Table 2 each have a repeated invariant
$\sigma$. As noted in Section 2.2, Conjecturing does not allow repeated
invariants and so these equations will not be recoverable as bounds. In
Section 5.4, we describe ways to address this deficiency and recover equations
such as I.6.2 and I.6.2.b.

The normalized root-mean-square error (NRMSE) calculated for 100 test examples
for the best-performing conjecture on the training data based on mean absolute
error. NRMSE is calculated as $1/\sigma_{f}$ multiplied by the root-mean-
square error on the test examples, where $\sigma_{f}$ is the standard
deviation of the invariant of interest for the test examples. Conjecturing
produced conjectures with NRMSE less than 1.0 for 7 of the 10 equations.

Despite the fact that Conjecturing is not designed for recovery of equations,
we see that it can be successful in doing so for lower-complexity, albeit
nonlinear, equations.

Table 6 in the Appendix contains the results of applying Conjecturing to the
datasets when noise is added to the invariant of interest. The noise in the
table is the standard deviation of the normal distribution with mean zero
added to the invariant of interest; the noise is the noise level tolerated by
AI Feynman as reported by Udrescu and Tegmark (2020). The time is the time
reported by Udrescu and Tegmark (2020) for AI Feynman to recover the equation.

Conjecturing is unable to achieve exact recovery of the equations with the
introduction of noise. The NRMSE is less than 1.0 for six of the 10 equations
and does not exceed 1.625. The strict requirements of the Dalmatian heuristic
seem to prevent exact recovery in the presence of noise, but Conjecturing is
still able to produce good approximations of the equations.

Performance of AI Feynman (Udrescu and Tegmark 2020) on Gravity and Real
Estate Examples. We apply the implementation of AI Feynman available here:
https://github.com/SJ001/AI-Feynman to the gravity example described in
Section 3.1 and the real estate example described in Section 3.2. A difference
between our gravity example and the datasets used by Udrescu and Tegmark
(2020) is that for the gravity example, the gravitational constant $G$ is the
same for every data point, but for the datasets used by Udrescu and Tegmark
(2020), $G$ is treated as a variable and is different for each point.

For the real estate data, we apply AI Feynman to the data to attempt to
discover the relationship between the property of interest and the input
features. The original numeric features are supplied along with boolean
features corresponding to the levels of the propertyType feature. Note that AI
Feynman is designed for recovering numeric functions and is therefore not
suitable for boolean relationships such as those in the real estate example.

For both instances, the AI Feynman implementation aborts with an error
regarding an eigenvalue calculation. We suspect that the source of the failure
in both cases may be due to the difference in treatment of constants. In our
gravity example, the gravitational constant is the same for all points while
in analogous examples, Udrescu and Tegmark (2020) treat constants as variables
and generate a unique value for each observation. Our practice of treating the
gravitational constant as the same for all observations may be contributing to
an error in matrix calculations for AI Feynman. In the real estate example,
each observation has a feature with the same value (the $300,000 cutoff). This
constant column in the data matrix could also be contributing to an error in
matrix calculations for AI Feynman. In the electronic companion, we include
the code and output for AI Feynman applied to 1) their Example 1,
demonstrating that our installation is functional, 2) our gravity example,
including the error message, and 3) our real estate example, including the
error message. These examples show that the Conjecturing framework can provide
useful insights on examples where AI Feynman cannot.

Table 2: Results for Conjecturing on Datasets from (Udrescu and Tegmark 2020) |  | Number of |  |  |  | Recovered
---|---|---|---|---|---|---
Instance | Equation | Invariants | Complexity | Time (s) | NRMSE |  by Eureqa?
I.6.2.a | $f=e^{-\theta^{2}/2}/\sqrt{2\pi}$ | 2 | 9 | $10^{3}$ | 0.000 | No
I.6.2 | $f=e^{-\theta^{2}/2\sigma^{2}}/\sqrt{2\pi\sigma^{2}}$ | 3 | 13 | NA | 0.553 | No
I.6.2.b | $f=e^{-(\theta-\theta_{1})^{2}/2\sigma^{2}}/\sqrt{2\pi\sigma^{2}}$ | 5 | 16 | NA | 1.092 | No
I.8.14 | $d=\sqrt{(x_{2}-x_{1})^{2}+(y_{2}-y_{1})^{2}}$ | 4 | 10 | $10^{3}$ | 0.000 | No
I.9.18 | $\frac{Gm_{1}m_{2}}{(x_{2}-x_{1})^{2}+(y_{2}-y_{1})^{2}+(z_{2}-z_{1})^{2}}$ | 9 | 17 | NA | 1.100 | No
I.10.7 | $m=\frac{m_{0}}{\sqrt{1-\frac{v^{2}}{c^{2}}}}$ | 3 | 9 | 5 | 0.000 | No
I.11.19 | $A=x_{1}y_{1}+x_{2}y_{2}+x_{3}y_{3}$ | 6 | 11 | NA | 1.014 | Yes
I.12.1 | $F=\mu N_{n}$ | 2 | 3 | 5 | 0.000 | Yes
I.12.2 | $F=\frac{q_{1}q_{2}}{4\pi\epsilon r^{2}}$ | 5 | 12 | NA | 0.671 | Yes
I.12.4 | $E_{f}=\frac{q_{1}}{4\pi\epsilon r^{2}}$ | 4 | 11 | NA | 0.135 | Yes

###  5.4 Experiments with the Nguyen Benchmark Suite Nguyen et al. (2011)

We apply our invariant conjecturing implementation to the Nguyen benchmark
suite (Nguyen et al. 2011) so as to draw comparisons with symbolic regression
methods described by Petersen et al. (2021). The Nguyen benchmark suite is a
set of 12 equations. As mentioned before, our Conjecturing framework is
designed for discovering nonlinear relationships in the form of bounds and
boolean relationships while symbolic regression methods are designed to
recover equations. In these experiments, we investigate the ability of
invariant conjecturing to recover equations (or approximations) among the
discovered bounds.

The benchmark equations are in Table 3. We generate an instance for each using
the protocols described by Petersen et al. (2021). For each equation, 20
training examples and 20 test examples are generated. For equations 1 though
6, $x$ is sampled from a Uniform(-1,1) distribution; for equation 7, $x$ is
sampled from a Uniform(0,2) distribution; for equation 8, $x$ is sampled from
a Uniform(0,4) distribution; for equations 9 through 12, $x$ and $y$ are
sampled from a Uniform(0,1) distribution. For each equation, we allow a time
limit of 10,000 seconds for generating upper and lower bounds. The operators
include unary operators sine, cosine, natural log, and natural exponential,
and binary operators addition, subtraction, multiplication, and division.

As noted in Section 2.2, our Conjecturing framework does not allow repeated
invariants in conjectures. Therefore, most of the equations in Table 3 cannot
be recovered by our framework. For each equation, we first evaluate the
ability of the conjectured bounds to approximate the equation by reporting the
normalized root mean squared error (NRMSE) as defined and reported by Petersen
et al. (2021) and described in Section 5.3. We report NRMSE for the conjecture
with the lowest mean average error for the training examples.

The results in Table 3 indicate that our Conjecturing framework is able to
recover only equation 11 $f=x^{y}$. It is able to do so despite the fact that
the exponent operator is not included. Conjecturing produces the expression
$e^{y\log(x)}$ and simplified it to $x^{y}$. The NRMSE values for the best
bounds for the other instances range from 0.22 to 1.08. These values are
larger than those reported for symbolic regression methods as reported in
Table 10 of Petersen et al. (2021); the methods include deep symbolic
regression (DSR) Petersen et al. (2021), priority queue training (PQT)
Abolafia et al. (2018), vanilla policy gradient (VPG) Petersen et al. (2021),
genetic programming (GP), and a method implemented in Mathematica based on
Markov chain Monte Carlo and nonlinear regression.

For expressions 1 through 8, the only invariants are the invariant of interest
$f$ and the input invariant $x$. Because of the fact that no invariants can be
repeated, the Conjecturing framework is limited to the application of only
unary operators to $x$ and no expressions with binary operators are produced.
As an example, for the first equation, the best conjectured lower bound is

$f(x)\geq\sin(e^{e^{\cos(\sin(\sin(\sin(\sin(\log(\sin(\sin(\cos(\cos(\sin(e^{x})))))))))))}}).$

We now describe how our Conjecturing framework can be adapted to allow for
repeated invariants, and report results for the adapted method. To address
repeated invariants, we can add additional invariants using commonly-occurring
functional forms. For equations 1 through 8, we add invariants $x^{2}$,
$x^{3}$, $x^{4}$, $x^{5}$, $x^{6}$, $\sin(x)$, $\cos(x)$, $\sqrt{x}$ along
with two copies of the constant 1, and the constant 2. Two copies of the
constant 1 are included because it appears twice in equation 7. For equations
9 through 12 we also add $y^{2}$ as an invariant. Recall that in our
implementation, while invariants cannot repeat in a conjectures, operators
can. Therefore, an alternative approach to addressing the constants is to
include the unary operator of addition by 1.

As shown in Table 3, the Conjecturing framework is able to exactly recover
each equation as a bound when the additional invariants are included so that
the NRMSE values are 0.000 for all equations.

The practice of adding the invariants that are nonlinear functions of the
original input might appear to be impractical. However, as suggested by
Brunton et al. (2016), specifying these invariants can reflect expert
knowledge on the system being investigated. They note that identifying
candidate functions for SINDy “must be a coordinated effort to incorporate
expert knowledge, feature extraction, and other advanced methods.”
Conjecturing offers distinct capabilities for discovery, as nonlinear
functions can be specified as invariants or may still be discovered so long as
they do not involve repeated input invariants.

Table 3: Results for Conjecturing on the Nguyen Benchmark Suite Nguyen et al. (2011) |  | Without Additional Invariants | With Additional Invariants
---|---|---|---
Instance | Equation | Recovered? | NRMSE | Recovered? | NRMSE
1 | $f=x^{3}+x^{2}+x$ | No | 0.94 | Yes | 0.00
2 | $f=x^{4}+x^{3}+x^{2}+x$ | No | 0.86 | Yes | 0.00
3 | $f=x^{5}+x^{4}+x^{3}+x^{2}+x$ | No | 1.01 | Yes | 0.00
4 | $f=x^{6}+x^{5}+x^{4}+x^{3}+x^{2}+x$ | No | 0.63 | Yes | 0.00
5 | $f=\sin{(x^{2})}\cos{(x)}-1$ | No | 0.22 | Yes | 0.00
6 | $f=\sin{(x)}+\sin{(x+x^{2})}$ | No | 0.44 | Yes | 0.00
7 | $f=\log{(x+1)}+\log{(x^{2}+1)}$ | No | 0.58 | Yes | 0.00
8 | $f=\sqrt{x}$ | No | 1.08 | Yes | 0.00
9 | $f=\sin{(x)}+\sin{(y^{2})}$ | No | 1.03 | Yes | 0.00
10 | $f=2\sin{(x)}\cos{(y)}$ | No | 0.60 | Yes | 0.00
11 | $f=x^{y}$ | Yes | 0.00 | Yes | 0.00
12 | $f=x^{4}-x^{3}+\frac{1}{2}y^{2}-y$ | No | 0.83 | Yes | 0.00

##  6 Application to COVID-19 Data

In this section, we demonstrate the Conjecturing framework on synthetic
patient-level COVID-19 data that was provided as part of the Veterans Health
Administration (VHA) Innovation Ecosystem and precisionFDA COVID-19 Risk
Factor Modeling Challenge (https://precision.fda.gov/challenges/11/view). The
data include synthetic veteran patient health records including medical
encounters, conditions, medications, and procedures. All subjects are located
in Massachusetts.

The goal of the challenge was to better understand risk and protective factors
for COVID-19 outcomes. Participants were asked to predict alive/deceased
status. Since our goal is to discover potential risk and protective factors,
we evaluate the performance of Conjecturing by checking the performance of the
feature relationships on holdout test data rather than on prediction accuracy.
Establishing the risk and protective factors as causal would require
additional controlled experiments.

Predictions were based on information obtained through December 31, 2019. In
the training data, we drop all information pertaining to events on or after
January 1, 2020 and drop subjects who died before January 1, 2020.

We add a number of features for each patient based on the data, described in
Table 7 in the Appendix. In addition, for each numeric observation, we created
invariants for the mean and most recent value and for each reported allergy,
device, immunization, procedure, and discretely-measured observation we create
a property corresponding to each level. In total, we use 309 invariants and
362 properties. We use a training set consisting of 100 subjects from each
outcome class (deceased/alive). We compare the results of applying
Conjecturing with classification and regression trees (CART) Breiman et al.
(1984) which is another interpretable method.

###  6.1 Results for Conjecturing

Upper and lower bounds are generated for each invariant, and for each outcome.
These bounds, along with the 362 properties in the data, are used as
properties for Conjecturing-PROP. Conjectures are generated for both outcomes.
The parameter $skips$ is set to 90%. We use the remaining 73,497 subjects as a
test data set thereby allowing us to ascertain the effects of potential
overfitting to the 200 subjects used for training.

Among those with COVID-19 in the test data, 5,468 (8.0%) have a status of
deceased, and 68,029 (92.0%) are alive. There are 38 conjectures for
sufficient conditions for alive status and 40 conjectures for sufficient
conditions for deceased status produced by the framework. Tables 8-11 in the
Appendix contain the conjectures and evaluations.

Tables 9 and 11 in the Appendix contain quantitative evaluations of the
performance of the conjectures. Each table contains the precision, support,
and lift of each conjecture. Note that each conjecture is a sufficient
condition expressed as a conditional statement. The precision is the
percentage of test examples for which the conditional statement evaluates to
true among those for which the antecedent is true. Precision may be thought of
as the “hit rate” of the conjecture. The support is the number of test
examples for which the antecedent evaluates to true. The lift is the ratio of
the precision to the proportion of examples for which the consequent is true.
If the lift is greater than 1, then the conjecture is better at identifying
people for which the consequent is true than a random selection from the
population.

Of the 38 conjectures for alive status, 22 (57.9%) have lift at least 1.00.
The lift ranges from 0.81 to 1.07; note that the maximum possible lift for a
conjecture for alive status is 1.08 (1/(68029/73497)). Of the 40 conjectures
for deceased status, 34 (85%) have lift at least 1.00. The lift ranges from
0.17 to 4.15.

Consider the sufficient conditions for deceased status in Table 10 in the
Appendix. The conjecture with the highest precision and lift is

$\textit{longitude}>-\textit{age}\times\textit{medicationsLifetimePercCovered}\rightarrow\textit{Deceased},$

and has a lift of 4.15 meaning that a subject for which
$\textit{longitude}>-\textit{age}\times\textit{medicationsLifetimePercCovered}$
is 4.15 times as likely to die as a randomly selected subject. The conjecture
indicates that subjects in the east who are older and have a larger percentage
of medications covered by the payer are at higher risk of death. The presence
of longitude in the conjecture could be an indication of higher risk in
population centers in the east such as Boston. The percentage of medications
covered by the VHA is higher for subjects with more preexisting conditions and
for those with more expensive medications because there is a low copay annual
cap (currently $700444https://www.va.gov/health-care/copay-rates/, accessed
July 10th, 2022). Further, the conjecture produces a suggestion of a
functional form for the relationship between these factors. The conjecture
confirms the CDC guidance that older subjects and those with more preexisting
conditions are a higher risk of death from
Covid555https://www.cdc.gov/coronavirus/2019-ncov/need-extra-
precautions/people-with-medical-conditions.html, accessed July 10th, 2022.

The conjecture with the second-highest precision and lift is

$\textit{medicationsActive}>\lfloor\textit{hemoglobinA1cHemoglobinTotalInBlood}\rfloor\rightarrow\textit{Deceased},$

with a lift of 3.29. The condition includes the number of active medications
and the ratio of hemoglobin A1c to total hemoglobin (an HbA1c test). The
conjecture indicates that those with more active medications than the HbA1c
percentage are at higher risk of death, which again agrees with the CDC
guidance concerning preexisting conditions. Typical values for HbA1c for non-
diabetic patients are below 5.7%, while diabetic subjects can have values
between 6.5% and 10.0% 666https://www.cdc.gov/diabetes/managing/managing-
blood-sugar/a1c.html, accessed July 10th, 2022. The number of active
medications can indicate a larger number of preexisting conditions, and the
conjecture suggests that for diabetic patients, those with additional
conditions are at higher risk.

The conjecture with the third-highest precision and lift is

$\textit{age}>\textit{carbonDioxide}\times\lfloor\textit{potassium}\rfloor\rightarrow\textit{Deceased},$

with a lift of 3.12. The conjecture suggests that older subjects with lower
$CO_{2}$ levels and lower potassium are at higher risk of death. Lower
$CO_{2}$ levels and abnormal potassium levels, particularly lower levels, has
been independently studied and associated with COVID-19 morbidity and
mortality (Hu et al. 2021, Noori et al. 2022). In addition to validating a
role for these invariants, the conjecture suggests a potential nonlinear
relationship among them.

For both outcomes, the Conjecturing framework is able to generate new
sufficient conditions that are true for the respective outcome at higher rates
than would be expected for a patient selected at random. These results
indicate that the conjecturing process is capturing relationships that hold
across the population and are not merely reflective of the 200 training
samples. In other words, overtraining appears to be mitigated. The discovered
relationships, and the direct and indirect relationships that they indicate
among features, are validated by the medical literature and provide
suggestions for deeper investigations into the functional form of the
relationships and the extent of causality. The number of conjectures generated
is not overwhelming for a human investigator to consider and further
investigate.

###  6.2 Comparison with an Interpretable Model

We now consider the results of applying classification and regression trees
(CART) (Breiman et al. 1984) to the COVID-19 data. A model is fit using the
implementation in the R library rpart (Therneau and Atkinson 2019). The tree
produced by rpart is

     n= 200

     node), split, n, loss, yval, (yprob)
           * denotes terminal node

      1) root 200 100 False (0.50000000 0.50000000)
        2) age< 66.48323 92  17 False (0.81521739 0.18478261)
          4) active_care_plan_length< 33.3039 68   5 False (0.92647059 0.07352941) *
          5) active_care_plan_length>=33.3039 24  12 False (0.50000000 0.50000000)
           10) Body_Height>=167.15 8   1 False (0.87500000 0.12500000) *
           11) Body_Height< 167.15 16   5 True (0.31250000 0.68750000) *
        3) age>=66.48323 108  25 True (0.23148148 0.76851852)
          6) mean_Potassium>=4.85 10   4 False (0.60000000 0.40000000) *
          7) mean_Potassium< 4.85 98  19 True (0.19387755 0.80612245) *

Nodes 4, 6, 7, 10, and 11, correspond to leaf nodes. Each leaf node
corresponds to a sufficient condition for deceased status (True) or alive
status (False). Nodes 4, 6, and 10 correspond to sufficient conditions for
alive status, and nodes 7 and 11 correspond to sufficient conditions for
deceased status. Tables 4 and 5 contain the conditions and quantitative
evaluations of the conditions produced by CART. When comparing the results
with those of Conjecturing in Tables 8-11 in the Appendix, we note that CART

  1. 1.

produces many fewer conditions (3 versus 38 for alive status, 2 versus 40 for
deceased status),

  2. 2.

produces conditions that are in conjunctive normal form where each clause
consists of a single numeric bound while Conjecturing tends to leverage
nonlinear relationships among invariants as the basis for conditions,

  3. 3.

produces two conditions with much larger support in the test data than those
produced by Conjecturing (node 4 has 36,531 and node 7 has 12,443),

  4. 4.

produces only two conditions that have lift greater than 1.0 (node 4 has lift
1.06 and node 7 has lift 2.85), and

  5. 5.

does not produce conditions with better precision or lift than the best
conditions produced by Conjecturing.

We note that both CART and Conjecturing are able to leverage categorical
variables for conditions, though CART does not do so for this training set. An
example of such a condition is conjecture 27 in Table 10 in the Appendix.

Similar to many decision tree frameworks, CART leverages univariate bounds as
component properties in its invariant clauses. Conjecturing is unlikely to
derive numeric bounds for individual features but instead produces more
nonlinear relationships between invariants. Decision tree frameworks such as
CART and Conjecturing are complementary approaches for discovery of patterns
among numeric and categorical features, but we see that Conjecturing is
capable of producing more complex yet interpretable relationships.

Table 4: Conditions from CART for Alive/Deceased Status Among Those with COVID $\begin{array}[]{lr}\textrm{Node Number}&\textrm{Sufficient Condition}\\\ \hline\cr\hline\cr 4&\textit{age}<66.48\ \&\ \textit{activeCarePlanLength}<33.30\rightarrow\textit{Alive}\\\ 6&\textit{age}\geq 66.48\ \&\ \textit{meanPotassium}\geq 4.85\rightarrow\textit{Alive}\\\ 10&\textit{age}<66.48\ \&\ \textit{activeCarePlanLength}\geq 33.30\ \&\ \textit{bodyHeight}\geq 167.15\rightarrow\textit{Alive}\\\ \hline\cr 7&\textit{age}\geq 66.48\ \&\ \textit{meanPotassium}<4.85\rightarrow\textit{Deceased}\\\ 11&\textit{age}<66.48\ \&\ \textit{activeCarePlanLength}\geq 33.30\ \&\ \textit{bodyHeight}<167.15\rightarrow\textit{Deceased}\\\ \hline\cr\end{array}$ Table 5: Evaluation of Conditions from CART Among Those with COVID Node Number | Consequent | Precision | Support | Lift
---|---|---|---|---
4 | Alive | 98.28% | 36531 | 1.06
6 | Alive | 81.69% | 2070 | 0.88
10 | Alive | 91.57% | 5371 | 0.99
7 | Deceased | 21.17% | 12443 | 2.85
11 | Deceased | 5.37% | 3746 | 0.72

##  7 Conclusions

We have demonstrated that automated search for conjectured feature-relations
can enhance existing learning algorithms. The discovery of these kinds of
feature relationships can also initiate new collaboration with domain
scientists and lead to new scientific knowledge.

Our Conjecturing framework was able to recover the functional form for gravity
with only the measured force, masses, and distance. The framework also
recovered a hidden interaction between price per square foot, square footage,
and price in real estate data that leads to improved classification
performance. Using synthetic patient-level COVID-19 data, the framework
produced conjectures that provide insight into the risk of death.

The current version of Conjecturing requires that conjectures are true for
every example. Future research will adjust the algorithm to better handle
noisy data by generating conjectures that do not necessarily hold for all
examples. If the handling of noise can be improved, then Conjecturing may be
able to be adapted to support predictive modeling efforts.

If the Conjecturing framework can provide functional relationships without
constants of proportionality, the constant can be determined using regression
with the original data. Suppose that the Conjecturing framework indicates a
relationship between the response $y$ and predictors $x$ of the form $y\leq
b_{1}f(x)$ for an unknown constant $b_{1}$. A regression model can be fit of
the form $\hat{y}=b_{0}+b_{1}f(x)$ using the data $(x_{i},y_{i})$,
$i=1,\ldots,n$. The best strategy for determining constants of proportionality
is another avenue for future research.

Another area for potential research involves the so-called $p>>n$ problem.
That is, if the number of features is larger than the number of observations,
then there are insufficient degrees of freedom to estimate a linear model with
all $p$ features or any more complex model. In such situations, feature and/or
model selection tools are needed to search over potentially large model
spaces. As desired model complexity increases (e.g. consideration of
interaction terms), searching over such large model spaces can become
computationally prohibitive. For instance, suppose an investigator seeks to
identify a model by selecting the “best” subset from among 10 features and
their associated 45 two-way interactions. In this example, simply considering
models with only 10 variables requires searching over a model space larger
than 2.9 billion. Future research will investigate the ability of the
Conjecturing framework to simplify model spaces and hence, provide a mechanism
for a more expeditious search of plausible models.

##  8 Additional Results for Conjecturing for Recovery of Equations

Table 6: Additional Results for Conjecturing on Datasets from (Udrescu and Tegmark 2020) |  | Number of |  |  |  |
---|---|---|---|---|---|---
Instance | Equation | Invariants | Complexity | Noise | Time (s) | NRMSE
I.6.2.a | $f=e^{-\theta^{2}/2}/\sqrt{2\pi}$ | 2 | 9 | $10^{-2}$ | 16 | 0.234
I.6.2 | $f=e^{-\theta^{2}/2\sigma^{2}}/\sqrt{2\pi\sigma^{2}}$ | 3 | 13 | $10^{-4}$ | 2992 | 0.873
I.6.2.b | $f=e^{-(\theta-\theta_{1})^{2}/2\sigma^{2}}/\sqrt{2\pi\sigma^{2}}$ | 5 | 16 | $10^{-4}$ | 4792 | 1.625
I.8.14 | $d=\sqrt{(x_{2}-x_{1})^{2}+(y_{2}-y_{1})^{2}}$ | 4 | 10 | $10^{-4}$ | 544 | 1.402
I.9.18 | $\frac{Gm_{1}m_{2}}{(x_{2}-x_{1})^{2}+(y_{2}-y_{1})^{2}+(z_{2}-z_{1})^{2}}$ | 9 | 17 | $10^{-5}$ | 5975 | 0.977
I.10.7 | $m=\frac{m_{0}}{\sqrt{1-\frac{v^{2}}{c^{2}}}}$ | 3 | 9 | $10^{-4}$ | 14 | 0.004
I.11.19 | $A=x_{1}y_{1}+x_{2}y_{2}+x_{3}y_{3}$ | 6 | 11 | $10^{-3}$ | 184 | 1.076
I.12.1 | $F=\mu N_{n}$ | 2 | 3 | $10^{-3}$ | 12 | 0.001
I.12.2 | $F=\frac{q_{1}q_{2}}{4\pi\epsilon r^{2}}$ | 5 | 12 | $10^{-2}$ | 17 | 0.746
I.12.4 | $E_{f}=\frac{q_{1}}{4\pi\epsilon r^{2}}$ | 4 | 11 | $10^{-2}$ | 12 | 1.040

##  9 Features Generated for COVID-19 Data

Table 7: Feature Definitions for COVID-19 Data Feature Name |  Definition
---|---
$healthcareExpenses$ |  The total lifetime cost of healthcare to the patient
$healthcareCoverage$ |  The total lifetime cost of healthcare services that were covered by payers
$latitude$ |  Latitude of patient’s home address
$longitude$ |  Longitude of patient’s home address
$age$ |  Current age of patient
$numAllergies$ |  Number of ongoing patient allergies
$activeCarePlans$ |  Number of current care plans
$lifetimeCarePlans$ |  Number of lifetime care plans
$activeCarePlanLength$ |  Length of time under current care plans
$lifetimeCarePlanLength$ |  Total lifetime length under care plans
$activeConditions$ |  Number of current health conditions
$lifetimeConditions$ |  Number of lifetime health conditions
$activeConditionLength$ |  Amount of time since current health condition(s) diagnosis
$lifetimeConditionLength$ |  Amount of time since first diagnosis of a health condition
$deviceLifetimeLength$ |  Total length of time using a medical device (e.g. pacemaker)
$encountersCount$ |  Total number of encounters with a healthcare professional
$encountersLifetimeTotalCost$ |  Total lifetime cost of healthcare encounters
$encountersLifetimeBaseCost$ |  Total lifetime cost of healthcare encounters, not including any line item costs related to medications, immunizations, procedures, or other services
$encountersLifetimePayerCoverage$ |  Total lifetime cost of healthcare encounters that were covered by payers
$encountersLifetimePercCovered$ |  Percentage of lifetime cost of healthcare encounters that were covered by payer
$imagingStudiesLifetime$ |  Number of lifetime imaging diagnostics (e.g. MRI) performed on patient
$immunizationsLifetime$ |  Number of lifetime immunizations received by patient
$immunizationsLifetimeCost$ |  Total lifetime cost of all immunizations received by patient
$medicationsLifetime$ |  Number of lifetime medications prescribed
$medicationsLifetimeCost$ |  Total lifetime cost of medications
$medicationsLifetimePercCovered$ |  Percentage of lifetime medication cost coverered by payer
$medicationsLifetimeLength$ |  Total lifetime length on prescribed medications
$medicationsLifetimeDispenses$ |  Total lifetime number of prescription dispenses
$medicationsActive$ |  Number of current prescriptions
$proceduresLifetime$ |  Number of lifetime medical procedures (e.g. surgery) performed on patient
$proceduresLifetimeCost$ |  Total lifetime cost of all medical procedures performed on patient

##  10 Results for Conjecturing Applied to COVID Data

Table 8: Conjectures for Alive Status Among Those with COVID
$\begin{array}[]{lr}&\textrm{Conjecture}\\\ \hline\cr\hline\cr
1&\textit{activeConditionLength}>\textit{age}^{2}/\textit{latitude}\rightarrow\textit{Alive}\\\
2&\textit{medicationsLifetime}<-\textit{immunizationsLifetimeCost}+2\times\textit{proceduresLifetime}\rightarrow\textit{Alive}\\\
3&\textit{medicationsActive}<\min\\{\textit{sodium},\lfloor\textit{QOLS}\rfloor\\}\rightarrow\textit{Alive}\\\
4&\textit{activeCarePlans}<e^{\textit{medicationsLifetimePercCovered}}-\textit{immunizationsLifetime}\rightarrow\textit{Alive}\\\
5&\textit{activeCarePlanLength}<10^{\textit{encountersLifetimePercCovered}}\times\textit{activeCarePlans}\rightarrow\textit{Alive}\\\
6&\textit{lifetimeConditionLength}>\sqrt{\textit{QALY}}^{\textit{activeConditions}}\rightarrow\textit{Alive}\\\
7&\textit{lifetimeCarePlans}>\textit{encountersCount}/2\rightarrow\textit{Alive}\\\
8&\textit{lifetimeCarePlans}>\min\\{\textit{triglycerides},\lfloor\textit{lifetimeConditions}\rfloor\\}\rightarrow\textit{Alive}\\\
9&\textit{encountersCount}>\textit{activeConditions}+\textit{medicationsLifetimeCost}+1\rightarrow\textit{Alive}\\\
10&\textit{diabetes }\&\textit{
latitude}<\sqrt{\textit{totalCholesterol}}+\textit{meanCarbonDioxide}\rightarrow\textit{Alive}\\\
11&\textit{activeCarePlans}>\textit{medicationsLifetime}^{\textit{medicationsLifetimeCost}}\rightarrow\textit{Alive}\\\
12&\textit{lifetimeCarePlans}>\sqrt{\textit{enountersCount}}+\textit{QOLS}\rightarrow\textit{Alive}\\\
13&\textit{age}<\textit{lifetimeConditions}\times\log(\textit{latitude})\rightarrow\textit{Alive}\\\
14&\textit{anemiaDisorder }\&\textit{
activeConditionLength}<\min\\{\textit{ureaNitrogen},|\textit{activeCarePlanLength}|\\}\rightarrow\textit{Alive}\\\
15&\textit{lifetimeCarePlans}>\max\\{\textit{DALY},\lfloor\textit{potassium}\rfloor\\}\rightarrow\textit{Alive}\\\
16&\textit{healthcareCoverage}>\lfloor\textit{DALY}\rfloor\times\textit{lifetimeCarePlanLength}\rightarrow\textit{Alive}\\\
17&\textit{lifetimeConditionLength}>(\textit{encountersLifetimeTotalCost}-1)/\textit{proceduresLifetime}\rightarrow\textit{Alive}\\\
18&\textit{healthcareCoverage}>\textit{lifetimeConditionLength}^{2}/\textit{imagingStudiesLifetime}\rightarrow\textit{Alive}\\\
19&\textit{immunizationsLifetimeCost}>(\textit{bodyHeight}-1)^{\textit{immunizationsLifetime}}\rightarrow\textit{Alive}\\\
20&\textit{numAllergies}<\textit{activeCarePlanLength}-\textit{age}+1\rightarrow\textit{Alive}\\\
21&\textit{activeCarePlans}<\lfloor\textit{microalbuminCreatineRatio}\rfloor-\textit{proceduresLifetimeCost}\rightarrow\textit{Alive}\\\
22&\textit{latitude}<\sqrt{\textit{encountersLifetimeTotalCost}}-\textit{medicationsLifetime}\rightarrow\textit{Alive}\\\
23&\textit{bodyMassIndex40}\rightarrow\textit{Alive}\\\
24&\textit{activeCarePlanLength}>\textit{age}+\textit{proceduresLifetime}-1\rightarrow\textit{Alive}\\\
25&\textit{medicationsLifetimeLength}<2\times\textit{activeCarePlan}\times\textit{deviceLifetimeLength}\rightarrow\textit{Alive}\\\
26&\textit{carbonDioxide}>\textit{respiratoryRate}\times\lfloor\textit{hemoglobinA1cHemoglobinTotalInBlood}\rfloor\rightarrow\textit{Alive}\\\
27&\textit{osteoporosisDisorder }\&\textit{
lifetimeCarePlanLength}<\min\\{\textit{painSeverity},2\times\textit{activeCarePlanLength}\\}\rightarrow\textit{Alive}\\\
28&\textit{healthcareCoverage}<\sqrt{\textit{encountersLifetimePayerCoverage}}\times\textit{medicationsLifetime}\rightarrow\textit{Alive}\\\
29&\textit{respiratoryRate}<\textit{painSeverity}+\lceil\textit{leukocytesVolumeInBlood}\rceil\rightarrow\textit{Alive}\\\
30&\textit{meanPainSeverity}>\max\\{\textit{proceduresLifetime},\lfloor\textit{hemoglobinA1cHemoglobinTotalInBlood}\rfloor\\}\rightarrow\textit{Alive}\\\
31&\textit{prediabetes }\&\textit{
meanDiastolicBloodPressure}>\lfloor\textit{carbonDioxide}\rfloor+\textit{meanHeartRate}\rightarrow\textit{Alive}\\\
32&\textit{latitude}<\textit{encountersCount}\times\lfloor\textit{QOLS}\rfloor\rightarrow\textit{Alive}\\\
33&\textit{painSeverity}<\lfloor\textit{meanPainSeverity}\rfloor-1\rightarrow\textit{Alive}\\\
34&\textit{lifetimeCarePlans}>\max\\{\textit{DALY},\lfloor\textit{potassium}\rfloor\\}\rightarrow\textit{Alive}\\\
35&\textit{latitude}<2\times\textit{DALY}\times\textit{encountersLifetimePercCovered}\rightarrow\textit{Alive}\\\
36&\textit{activeCarePlanLength}>\max\\{\textit{sodium},\lceil\textit{hematocritVolume}\rceil\\}\rightarrow\textit{Alive}\\\
37&\textit{medicationsLifetimePercCovered}>\textit{latitude}^{2}/\textit{medicationsLifetimeDispense}\rightarrow\textit{Alive}\\\
38&\textit{healthcareCoverage}>\frac{\textit{healthcareExpenses}}{2\times\textit{encountersLifetimePercCovered}}\rightarrow\textit{Alive}\\\
\hline\cr\end{array}$

Table 9: Evaluation of Conjectures for Alive Status Among Those with COVID Conjecture | Precision | Support | Lift
---|---|---|---
1 | 99.11% | 3042 | 1.07
2 | 98.54% | 3700 | 1.06
3 | 98.51% | 9480 | 1.06
4 | 98.25% | 8231 | 1.06
5 | 98.23% | 6103 | 1.06
6 | 98.05% | 7937 | 1.06
7 | 97.79% | 4161 | 1.06
8 | 97.65% | 5453 | 1.05
9 | 97.65% | 2939 | 1.05
10 | 97.41% | 1969 | 1.05
11 | 97.07% | 4774 | 1.05
12 | 96.79% | 2242 | 1.05
13 | 96.61% | 1650 | 1.04
14 | 95.65% | 2045 | 1.03
15 | 95.39% | 6311 | 1.03
16 | 95.31% | 1236 | 1.03
17 | 95.20% | 1687 | 1.03
18 | 95.07% | 4017 | 1.03
19 | 94.32% | 440 | 1.02
20 | 94.00% | 500 | 1.02
21 | 93.59% | 1498 | 1.01
22 | 93.25% | 1185 | 1.01
23 | 92.09% | 834 | 0.99
24 | 91.01% | 1213 | 0.98
25 | 90.79% | 999 | 0.98
26 | 89.77% | 831 | 0.97
27 | 89.41% | 727 | 0.97
28 | 88.61% | 2389 | 0.96
29 | 88.27% | 358 | 0.95
30 | 88.21% | 704 | 0.95
31 | 87.67% | 1890 | 0.95
32 | 87.20% | 2250 | 0.94
33 | 87.14% | 583 | 0.94
34 | 85.79% | 1612 | 0.93
35 | 84.50% | 755 | 0.91
36 | 82.76% | 586 | 0.89
37 | 76.22% | 677 | 0.82
38 | 74.75% | 99 | 0.81

Table 10: Conjectures for Deceased Status Among Those with COVID
$\begin{array}[]{lr}\textrm{Conjecture}&\textrm{Sufficient Condition}\\\
\hline\cr\hline\cr
1&\textit{longitude}>-\textit{age}\times\textit{medicationsLifetimePercCovered}\rightarrow\textit{Deceased}\\\
2&\textit{medicationsActive}>\lfloor\textit{hemoglobinA1cHemoglobinTotalInBlood}\rfloor\rightarrow\textit{Deceased}\\\
3&\textit{age}>\textit{carbonDioxide}\times\lfloor\textit{potassium}\rfloor\rightarrow\textit{Deceased}\\\
4&\textit{deviceLifetimeLength}\leq
2\times\textit{creatinine}^{\textit{healthcareExpenses}}\rightarrow\textit{Deceased}\\\
5&\textit{implantableCardiacPacem}\rightarrow\textit{Deceased}\\\
6&\textit{latitude}<\log(\textit{age})/\log(10)^{\textit{activeCarePlans}}\rightarrow\textit{Deceased}\\\
7&\textit{medicationsActive}>\lceil\log(\textit{alkalinePhosphataseEnzymatic
Activity})/\log(10)\rceil\rightarrow\textit{Deceased}\\\
8&\textit{immunizationsLifetimeCost}<\textit{age}\times\textit{immunizationsLifetime}^{2}\rightarrow\textit{Deceased}\\\
9&\textit{colonoscopy }\&\textit{
coronaryHeartDisease}\rightarrow\textit{Deceased}\\\
10&\textit{activeCarePlans}<\min\\{\textit{deviceLifetimeLength},\textit{medicationsActive}\\}\rightarrow\textit{Deceased}\\\
11&\textit{glucose}>\lceil\textit{creatinine}\times\textit{meanGlucose}\rightarrow\textit{Deceased}\\\
12&\textit{bodyWeight}>\lfloor\textit{meanBodyWeight}\rfloor+1\rightarrow\textit{Deceased}\\\
13&\textit{healthcareExpenses}<\textit{deviceLifetimeLength}^{2}\times\textit{lifetimeConditionLength}\rightarrow\textit{Deceased}\\\
14&\textit{lifetimeConditions}>\textit{activeCarePlans}+\lfloor\textit{ureaNitrogen}\rfloor\rightarrow\textit{Deceased}\\\
15&\textit{activeCarePlans}>\lfloor\log(\textit{triglycerides})\rfloor\rightarrow\textit{Deceased}\\\
16&\textit{healthcareExpenses}<\textit{lifetimeConditionLength}^{2}+\textit{encountersLifetimeTotalCost}\rightarrow\textit{Deceased}\\\
17&\textit{overlappingMalignantNeo}\rightarrow\textit{Deceased}\\\
18&\textit{latitude}>\textit{ureaNitrogen}\lceil\textit{albuminMassVolumeInSerumOrPlasma}\rceil\rightarrow\textit{Deceased}\\\
19&\textit{activeConditionLength}>\textit{erythrocytesVolumeInBlood}\times\lceil\textit{hemoglobinMassVolumeInBlood}\rceil\rightarrow\textit{Deceased}\\\
20&\textit{chronicObstructiveBronc}\rightarrow\textit{Deceased}\\\
21&\textit{longitude}>\sqrt{\textit{healthcareCoverage}}-\textit{encountersCount}\rightarrow\textit{Deceased}\\\
22&\textit{age}>10^{\textit{medicationsActive}}-\textit{longitude}\rightarrow\textit{Deceased}\\\
23&\textit{chloride}<\lfloor\textit{meanChloride}\rfloor-\textit{lifetimeCarePlans}\rightarrow\textit{Deceased}\\\
24&\textit{medicationsActive}>\max\\{\textit{respiratoryRate},\log(\textit{latitude})\\}\rightarrow\textit{Deceased}\\\
25&\textit{localizedPrimaryOsteoa}\rightarrow\textit{Deceased}\\\
26&\textit{rheumatoidArthritis}\rightarrow\textit{Deceased}\\\
27&\textit{chronicPain }\&\textit{
smokesTobaccoDaily}\rightarrow\textit{Deceased}\\\
28&\textit{latitude}<\lfloor\textit{erythrocyteDistributionWidth}\rfloor-\textit{meanPainSeverity}\rightarrow\textit{Deceased}\\\
29&\textit{activeCarePlans}>10^{\textit{medicationsActive}}/\textit{imagingStudiesLifetime}\rightarrow\textit{Deceased}\\\
30&\textit{tubalPregnancy}\rightarrow\textit{Deceased}\\\
31&\textit{activeConditions}<\textit{medicationsActive}^{2}-\textit{medicationsLifetime}\rightarrow\textit{Deceased}\\\
32&\textit{alcoholism }\&\textit{
majorDepressionDisorder}\rightarrow\textit{Deceased}\\\
33&\textit{creatinine}<\lceil\textit{meanCreatinine}\rceil/\textit{lifetimeCarePlans}\rightarrow\textit{Deceased}\\\
34&\textit{healthcareCoverage}<\textit{encoutnersLifetimePayerCoverage}\times\log(\textit{latitude})/\log(10)\rightarrow\textit{Deceased}\\\
35&\textit{encountersCount}<\min\\{\textit{DALY},10^{\textit{immunizationsLifetime}}\\}\rightarrow\textit{Deceased}\\\
36&\textit{lifetimeCarePlanLength}>\textit{age}+e^{\textit{medicationsLifetime}}\rightarrow\textit{Deceased}\\\
37&\textit{healthcareCoverage}<\textit{activeCondionLength}^{2}-\textit{encountersLifetimeTotalCost}\rightarrow\textit{Deceased}\\\
38&\textit{age}>1/2\times\textit{healthcareExpenses}/\textit{immunizationsLifetimeCost}\rightarrow\textit{Deceased}\\\
39&\textit{activeCarePlanLength}>\textit{activeConditionLength}\times
e^{\textit{DALY}}\rightarrow\textit{Deceased}\\\
40&\textit{medicationsLifetime}<\sqrt{\textit{encountersLifetimePayerCoverage}}-\textit{age}\rightarrow\textit{Deceased}\\\
\hline\cr\end{array}$

Table 11: Evaluation of Conjectures for Deceased Status Among Those with COVID Conjecture | Precision | Support | Lift
---|---|---|---
1 | 30.91% | 372 | 4.15
2 | 24.44% | 2954 | 3.29
3 | 23.22% | 1722 | 3.12
4 | 22.47% | 632 | 3.02
5 | 22.39% | 844 | 3.01
6 | 21.74% | 1490 | 2.92
7 | 21.44% | 4427 | 2.88
8 | 21.24% | 2199 | 2.85
9 | 21.01% | 257 | 2.82
10 | 20.90% | 799 | 2.81
11 | 20.03% | 1058 | 2.69
12 | 18.25% | 548 | 2.45
13 | 17.56% | 467 | 2.36
14 | 17.54% | 1898 | 2.36
15 | 17.31% | 3328 | 2.33
16 | 17.23% | 940 | 2.32
17 | 16.92% | 130 | 2.27
18 | 16.13% | 1091 | 2.17
19 | 15.68% | 797 | 2.11
20 | 14.11% | 900 | 1.90
21 | 14.08% | 781 | 1.89
22 | 14.08% | 2230 | 1.89
23 | 11.99% | 884 | 1.61
24 | 11.84% | 1884 | 1.59
25 | 11.76% | 2159 | 1.58
26 | 11.44% | 201 | 1.54
27 | 10.54% | 607 | 1.42
28 | 10.27% | 1724 | 1.38
29 | 9.39% | 213 | 1.26
30 | 8.85% | 2147 | 1.19
31 | 8.22% | 4150 | 1.10
32 | 7.79% | 1129 | 1.05
33 | 7.77% | 2408 | 1.04
34 | 7.73% | 634 | 1.04
35 | 5.20% | 3209 | 0.70
36 | 3.89% | 4602 | 0.52
37 | 3.85% | 467 | 0.52
38 | 2.65% | 339 | 0.36
39 | 2.29% | 3188 | 0.31
40 | 1.30% | 1001 | 0.17

High Performance Computing resources provided by the High Performance Research
Computing (HPRC) Core Facility at Virginia Commonwealth University
(https://chipc.vcu.edu) were used for conducting the research reported in this
work.

## References

  * Abolafia et al. (2018) Abolafia D, Norouzi M, Shen J, Zhao R, Le Q (2018) Neural Program Synthesis with Priority Queue Training, http://arxiv.org, arXiv:1801.03526v2.pdf.
  * Aghaei et al. (2021) Aghaei S, Gómez A, Vayanos P (2021) Strong Optimal Classification Trees, http://arxiv.org, arXiv:2103.15965.pdf.
  * Akyüz and Birbil (2021) Akyüz MH, Birbil ŞI (2021) Discovering Classification Rules for Interpretable Learning with Linear Programming, http://arxiv.org, arXiv:2104.10751.pdf.
  * Bellomarini et al. (2020) Bellomarini L, Benedetto D, Gottlob G, Sallinger E (2020) Vadalog: A Modern Architecture for Automated Reasoning with Large Knowledge Graphs. _Information Systems_ in press, 101528.
  * Bertsimas and Dunn (2017) Bertsimas D, Dunn J (2017) Optimal Classification Trees. _Machine Learning_ 106:1039–1082.
  * Blanquero et al. (2021) Blanquero R, Carrizosa E, Molero-Río C, Romero Morales D (2021) Optimal Randomized Classification Trees. _Computers and Operations Research_ 132, 105281.
  * Bradford et al. (2020) Bradford A, Day J, Hutchinson L, Larson CE, Mills M, Muncy D, Kaperick B, Van Cleemput N (2020) Automated conjecturing II: Chomp and Intelligent Game Play. _Journal of Artificial Intelligence Research_ 68:447–461.
  * Breiman et al. (1984) Breiman L, Friedman J, Olshen R, CJ S (1984) _Classification and Regression Trees_ (Wadsworth & Brooks).
  * Brunton et al. (2016) Brunton S, Proctor J, Kutz J (2016) Discovering Governing Equations from Data by Sparse Identification of Nonlinear Dynamical Systems. _Proceedings of the National Academy of the Sciences_ 113:3932–3937.
  * Chattopadhyay and Lipson (2014) Chattopadhyay I, Lipson H (2014) Data Smashing: Uncovering Lurking Order in Data. _Journal of the Royal Society Interface_ 11:20140826.
  * Chvátal (1972) Chvátal V (1972) On Hamilton’s Ideals. _Journal of Combinatorial Theory, Series B_ 12:163–168.
  * Chvátal and Erdös (1972) Chvátal V, Erdös P (1972) A note on hamiltonian circuits. _Discrete Mathematics_ 2(2):111–113.
  * Dash et al. (2018) Dash S, Günlük O, Wei D (2018) Boolean Decision Rules via Column Generation. _The Thirty-Second Conference on Neural Information Processing Systems (NeurIPS-18)_.
  * Elton (2020) Elton D (2020) Self-Explaining AI as an Alternative to Interpretable AI, http://arxiv.org, arXiv:2002.05149v4.pdf.
  * Fajtlowicz (1995) Fajtlowicz S (1995) On Conjectures of Graffiti. V. _Graph Theory, Combinatorics, and Algorithms, Vol. 1_ , 367–376, Wiley-Intersci. Publ. (New York: Wiley).
  * Fürnkranz et al. (2020) Fürnkranz J, Kliegr T, Paulheim H (2020) On Cognitive Preferences and the Plausibility of Rule-Based Models. _Machine Learning_ 109:853–898.
  * Haemers (1979) Haemers W (1979) On Some Problems of Lovász Concerning the Shannon Capacity of a Graph. _IEEE Transactions on Information Theory_ 25(2):231–232.
  * Hammer and Bonates (2006) Hammer P, Bonates T (2006) Logical Analysis of Data - An Overview: From Combinatorial Optimization to Medical Applications. _Annals of Operations Research_ 148\.
  * Hu et al. (2021) Hu D, Li J, Gao R, Wang S, Li Q, Chen S, Huang J, Huang Y, Li M, Long W, Liu Z, Guo L, Wu X (2021) Decreased CO2 Levels as Indicators of Possible Mechanical Ventilation-Induced Hyperventilation in COVID-19 Patients: A Retrospective Analysis. _Frontiers in Public Health_ 8:596168.
  * Jantzen (2020) Jantzen B (2020) Dynamical Kinds and Their Discovery, http://arxiv.org, arXiv:1612.04933.pdf.
  * Kanter and Veeramachaneni (2015) Kanter JM, Veeramachaneni K (2015) Deep Feature Synthesis: Towards Automating Data Science Endeavors. _2015 IEEE International Conference on Data Science and Advanced Analytics_.
  * Katz et al. (2016) Katz G, Shin ECR, Song D (2016) ExploreKit: Automatic Feature Generation and Selection. _16th IEEE International Conference on Data Mining_.
  * Khurana et al. (2018) Khurana U, Samulowitz H, Turaga D (2018) Feature Engineering for Predictive Modeling Using Reinforcement Learning. _The Thirty-Second AAAI Conference on Artificial Intelligence (AAAI-18)_.
  * Langely et al. (1987) Langely P, Simon HA, Bradshaw GL, Zytkow JM (1987) _Scientific Discovery: Computational Explorations of the Creative Process_ (MIT Press, Cambridge, MA).
  * Langley (2019) Langley P (2019) Scientific Discovery, Causal Explanation, and Process Model Induction. _Mind & Society_ 18:43–56.
  * Larson and Van Cleemput (2016) Larson CE, Van Cleemput N (2016) Automated conjecturing I: Fajtlowicz’s Dalmatian heuristic revisited. _Artificial Intelligence_ 231:17–38.
  * Larson and Van Cleemput (2017) Larson CE, Van Cleemput N (2017) Automated conjecturing III: Property-relations conjectures. _Annals of Mathematics and Artificial Intelligence_ 81(3):315–327.
  * Lovász (1979) Lovász L (1979) On the Shannon capacity of a graph. _Information Theory, IEEE Transactions on_ 25(1):1–7.
  * Lu et al. (2019) Lu J, Lee DK, Kim T, Danks D (2019) Good Explanation for Algorithmic Transparency. _Available at SSRN: http://dx.doi.org/10.2139/ssrn.3503603_ .
  * Nguyen et al. (2011) Nguyen Q, Nguyen X, O’Neill M, McKay R, Galván-López (2011) Semantically-Based Crossover in Genetic Programming: Application to Real-Valued Symbolic Regression. _Genetic Programming and Evolvable Machines_ 12:91–119.
  * Nicolau and Agapitos (2021) Nicolau M, Agapitos A (2021) Choosing Function Sets with Better Generalisation Performance for Symbolic Regression Models. _Genetic Programming and Evolvable Machines_ 22:73–100.
  * Noori et al. (2022) Noori M, Nejadghaderi S, Sullman M, Carson-Chahhoud K, Kolahi AA, Safiri S (2022) Epidemiology, Prognosis and Management of Potassium Disorders in Covid-19. _Reviews in Medical Virology_ 32:e2262.
  * Petersen et al. (2021) Petersen B, Larma M, Mundhenk T, Santiago C, Kim S, Kim J (2021) Deep Symbolic Regression: Recovering Mathematical Expressions from Data via Risk-Seeking Policy Gradients. _Proceedings of the International Conference on Learning Representation (ICLR)_.
  * Rudin (2019) Rudin C (2019) Stop Explaining Black Box Machine Learning Models for High Stakes Decisions and Use Interpretable Models Instead. _Nature Machine Intelligence_ 1:206–215.
  * Rudin and Ertekin (2018) Rudin C, Ertekin S (2018) Learning Customized and Optimized Lists of Rules with Mathematical Programming. _Mathematical Programming Computation_ 10:659–702.
  * Samek and Müller (2019) Samek W, Müller KR (2019) Towards Explainable Artificial Intelligence. _Explainable AI: Interpreting, Explaining and Visualizing Deep Learning_ , 5–22, Samek, W. and Montavon, G. and Vedaldi, A. and Hanson, L. and Müller, K. R. (eds) (Springer).
  * Schmidt and Lipson (2009) Schmidt M, Lipson H (2009) Distilling Free-Form Natural Laws from Experimental Data. _Science_ 324(5923):81–85.
  * Schrijver (2003) Schrijver A (2003) _Combinatorial optimization: polyhedra and efficiency_ , volume 24 (Springer Science & Business Media).
  * Tallorin et al. (2018) Tallorin L, Wang JL, Kim WE, Sahu S, Kosa NM, Yang P, Thompson M, Gilson MK, Frazier PI, Burkart MD, Gianneschi NC (2018) Discovering de novo Peptide Substrates for Enzymes using Machine Learning. _Nature Communications_ 9(1):1–10.
  * Therneau and Atkinson (2019) Therneau T, Atkinson B (2019) _rpart: Recursive Partitioning and Regression Trees_. URL https://CRAN.R-project.org/package=rpart, R package version 4.1-15.
  * Tibshirani (1996) Tibshirani R (1996) Regression Shrinkage and Selection via the LASSO. _Journal of the Royal Statistical Society: Series B_ 58:267–288.
  * Tsang et al. (2018a) Tsang M, Cheng D, Liu Y (2018a) Detecting Statistical Interactions from Neural Network Weights. _The Sixth International Conference on Learning Representations (ICLR-18)_.
  * Tsang et al. (2018b) Tsang M, Liu H, Purushotham S, Pavankumar M, Liu Y (2018b) Neural Interaction Transparency (NIT): Disentangling Learned Interactions for Improved Interpretability. _The Thirty-Second Conference on Neural Information Processing Systems (NeurIPS-18)_.
  * Tsang et al. (2020) Tsang M, Rambhatla H, Liu Y (2020) How Does This Interaction Affect Me? Interpretable Attribution for Feature Interactions. _The Thirty-Fourth Conference on Neural Information Processing Systems (NeurIPS-20)_.
  * Udrescu and Tegmark (2020) Udrescu SM, Tegmark M (2020) AI Feynman: A Physics-Inspired Method for Symbolic Regression. _Science Advances_ 6, eaay2631.
  * Verwer and Zhang (2019) Verwer S, Zhang Y (2019) Learning Optimal Classification Trees Using a Binary Linear Program Formulation. _The Thirty-Third AAAI Conference on Artificial Intelligence (AAAI-19)_.
  * Vilone and Longo (2020) Vilone G, Longo L (2020) Explainable Artificial Intelligence: A Systematic Review, http://arxiv.org, arXiv:2006.00093.pdf.
  * Wang and Rudin (2015) Wang F, Rudin C (2015) Falling Rule Lists. _The 18th International Conference on Artificial Intelligence and Statistics (AISTATS)_.
  * Wang et al. (2017) Wang T, Rudin C, Doshi-Velez F, Liu Y, Klampfl E, MacNeille P (2017) A Bayesian Framework for Learning Rule Sets for Interpretable Classification. _Journal of Machine Learning Research_ 18:1–37.
  * West (2001) West DB (2001) _Introduction to Graph Theory, $2^{nd}$ edition_ (Prentice Hall).