# Big City Bias: Evaluating the Impact of Metropolitan Size on Computational
Job Market Abilities of Language Models

Charlie Campanella
Indeed
<EMAIL_ADDRESS>
Rob van der Goot
IT University of Copenhagen
<EMAIL_ADDRESS>

###### Abstract

Large language models (LLMs) have emerged as a useful technology for job
matching, for both candidates and employers. Job matching is often based on a
particular geographic location, such as a city or region. However, LLMs have
known biases, commonly derived from their training data. In this work, we aim
to quantify the metropolitan size bias encoded within large language models,
evaluating zero-shot salary, employer presence, and commute duration
predictions in 384 of the United States’ metropolitan regions. Across all
benchmarks, we observe negative correlations between the metropolitan size and
the performance of the LLMS, indicating that smaller regions are indeed
underrepresented. More concretely, the smallest 10 metropolitan regions show
upwards of 300% worse benchmark performance than the largest
10.111https://github.com/charlie-campanella/big-city-bias

Big City Bias: Evaluating the Impact of Metropolitan Size on Computational Job
Market Abilities of Language Models

Charlie Campanella Indeed<EMAIL_ADDRESS>Rob van der Goot IT
University of Copenhagen<EMAIL_ADDRESS>

##  1 Introduction

Recent large language models (LLMs) are primarily trained on internet-derived
corpora Brown et al. (2020); Touvron et al. (2023a). These underlying datasets
are prone to linguistic and geographic bias. For example, the training data of
Common Crawl, used to train OpenAI’s GPT-3 Brown et al. (2020) and Meta AI’s
Llama Touvron et al. (2023a), is composed of 46% English language documents.
Such lexical imbalances contribute to an anglophone bias in various tasks,
exemplified by GPT-3.5’s "English-first" approach when translating
Irish/Gaeilge Chiarain et al. (2023) and its inability to pass Indonesian
primary school exams Koto et al. (2023) while simultaneously passing U.S.
college-entry exams OpenAI et al. (2023). Evidently, under-representation in
training corpora can adversely affect language model performance across
various tasks and contexts.

A significant population disparity exists among metropolitan regions in the
United States, with the largest, New York–Newark–Jersey City, NY-NJ, having
over 300 times more residents than the smallest, Eagle Pass, TX. These
population disparities seem to correlate with the amount of associated
information available online. For instance, querying “New York, NY” on
Wikipedia yields 114,067 results while a search for “Eagle Pass, TX” only
returns 1,516 results. The same queries on Google result in 1.7 billion and
8.2 million results, respectively. This training data disparity prompts us to
consider whether language models exhibit comparable biases in performance,
excelling in tasks associated with larger metropolitan regions.

In this work, we quantify the "big city bias" by evaluating salary, employer
presence, and commute duration predictions in the 384 metropolitan statistical
areas (MSAs) defined by the United States Census Bureau. Across all
benchmarks, analysis indicates a correlation between MSA population and
predictive accuracy, indicating superior language model efficacy in the
context of larger cities. Given the concentration of technical talent in the
United States’ 10 most populous metropolitan areas
222https://www.bls.gov/oes/current/oes151252.htm, these findings might inspire
AI practitioners to look beyond so-called "tech hubs" when applying LLM-
generated synthetic data to geographically diverse job matching applications.

We contribute:

  * •

A dataset with population, employer presence, and commute duration for 384
metropolitan areas in the USA.

  * •

Outputs of 5 recent language models predicting employer presence, commute
duration, and salary.

  * •

An analysis of the results of the language models for this data and task,
showing that larger metropolitan areas are more closely associated with
accurate predictive outcomes.

##  2 Methodology

###  2.1 Data and tasks

Job matching systems are designed to find the best-fitting role for a given
candidate. For example, Indeed.com, the world’s most popular job site
333https://www.indeed.com/about, considers a role to be a "bad fit" if there
is a salary or geographic mismatch,
444https://engineering.indeedblog.com/blog/2019/09/jobs-filter/ among other
criterion. Our geographically-focused evaluation prompts language models to
make salary, employer presence, and commute duration predictions for each
metropolitan region to evaluate their performance at common matching tasks.

Prediction biases can reduce the effectiveness of LLM-powered job matching
systems. For example, a system which fails to predict salaries within a given
region may match job seekers with financially incompatible roles. Inaccurate
commute duration forecasts can lead job seekers to consider positions that are
logistically unfeasible. Furthermore, not discerning whether an employer’s
regional presence is a major facility or a branch office can mislead job
seekers about career progression prospects.

As mentioned in Section 1, this analysis is conducted on metropolitan areas as
defined by the United States Census Bureau. For each task, we leverage data
collected from a reputable source, with distributions for each target variable
plotted in Figure 1. The data for each task is described in detail below.

Figure 1: For each variable the cumulative number of metropolitan areas which
match a certain threshold (y-axis). Note that the target variables (all except
population) are an average over multiple instances within the region, and are
thus not directly comparable (i.e. they could be over different sets of jobs,
employers, or commutes).

#### Salary Task Data:

Sourced from a confidential and proprietary dataset by Indeed.com which
includes pre-tax, base annual salaries for registered nurses, software
engineers, and warehouse workers in each MSA.

#### Employer Presence Task Data:

Based on a publicly accessible dataset published by People Data Labs
555https://www.peopledatalabs.com/top-employers-dataset which presents the top
10 employers and their respective employment counts across 384 U.S.
metropolitan regions. Our modified version of this dataset incorporates
additional "state_code" and "metro_state_code" columns, aligning with MSA
naming conventions used by the U.S. Census Bureau for standardized
referencing.

#### Commute Duration Task Data:

Derived from a custom script which uses the Google Maps API to generate 5
random commute origins, destinations, and durations within 20km of the
geographic center for each MSA. Note that each duration estimate assumes
driving as the mode of transport.

We plot the correlation of the population size and each of the target
variables against each other in Figure 2. From this, we can clearly see that
the average employer presence and MSA population have a very strong
correlation; this makes intuitive sense, as there is greater possibility in
larger cities to have larger employers. Other combinations of variables,
especially salary and commute duration, show less significant correlation,
although the slopes in Figure 1 have similar shapes.

Figure 2: Pearson correlations between MSA population and each target
variable. Note that the bottom-left and top-right are mirrored as Pearson
correlations are not directional.

###  2.2 Models

Our selection of language models was optimized for diversity in architectures,
training methods, and parameter size while minimizing computational costs.
Based on these criterion, we selected mistral-instruct-7b-v0.1 Jiang et al.
(2023), llama-2-chat-7b, llama-2-chat-70b Touvron et al. (2023b),
gpt-3.5-turbo Brown et al. (2020), and gpt-4 OpenAI et al. (2023).

We employ these models using prompting, designing task-specific prompts which
contain: the question, relevant entities for the task, and instructions on how
to structure the answer, so that the numerical response can be automatically
extracted. All prompts are shown below:

Prompt 2.1: Salary Predictions What is the average annual salary for a
{JOB_TITLE} in {METRO_AREA}? Return one estimate with no other text, do not
return a range. Salary:

Prompt 2.2: Employer Presence Predictions How many people does {EMPLOYER_NAME}
employ in {METRO_AREA}? Provide an estimate even if this information is not
publicly available. Return one estimate with no other text, do not return a
range. Number employed:

Prompt 2.3: Commute Duration Predictions What is the average driving commute
time from {ORIGIN} to {DESTINATION} in {METRO_AREA}? Provide an estimate even
though you are not a GPS. Return an estimate in minutes with no other text.
Commute Estimate:

###  2.3 Metrics

After obtaining a response from the language models (Section 2.2), we first
extract numerical values through the application of a regular expression
(`[^0-9\.]+$`). We compare these numeric values to the gold labels (Section
2.1), and then compute the percentage error for each task. For each task
category in a metropolitan region, we take into account 3-5 model outputs and
derive an average percentage error. This smooths out individual output
variances, ensuring a more reliable evaluation.

Next, we look at the Pearson correlation Pearson (1901) between the log of the
population size and the average prediction error. This will tell us whether
larger metropolitan areas consistently have more accurate predictions or vice
versa. We then employ the median error, as a quantitative metric of the
performance of the models. Since this is taken over the values of the
variables we predict, and they vary across tasks (Section 2.1), we use our
third metric: coefficient of determination (r2). This metric explains the
proportion of the variation in the error that is predictable from the
metropolitan size, and normally ranges from 0-1.

##  3 Results

LLM | Size | Correl. | Median | r2
---|---|---|---|---
|  |  | Err.% |
Salary |  |  |  |
mistral-instruct | 7b | -.0890∗ | 21.95 | .0079
llama-2-chat | 7b | -.3553∗ | 22.92 | .1262
llama-2-chat | 70b | -.2765∗ | 19.47 | .0765
gpt-3.5-turbo | 20b | -.3398∗ | 19.49 | .1155
gpt-4 | ? | -.4113∗ | 17.75 | .1692
Employer Presence |  |  |  |
mistral-instruct | 7b | -.2125∗ | 182.60 | .0452
llama-2-chat | 7b | -.3307∗ | 201.56 | .1094
llama-2-chat | 70b | -.2731∗ | 149.36 | .0746
gpt-3.5-turbo | 20b | -.3004∗ | 161.33 | .0902
gpt-4 | ? | -.3179∗ | 143.45 | .1010
Commute Duration |  |  |  |
mistral-instruct | 7b | -.1422∗ | 140.22 | .0202
llama-2-chat | 7b | -.1141∗ | 141.12 | .0130
llama-2-chat | 70b | -.1915∗ | 126.86 | .0367
gpt-3.5-turbo | 20b | -.0083∗ | 71.11 | .0001
gpt-4 | ? | -.1872∗ | 51.72 | .0351
Table 1: Pearson correlations between the log of the population and the average prediction error, median error, and coefficient of determination (r2). ∗ p $<$ 0.0001. Figure 3: The average error plotted against the log of the population for each task (shown on the right) and each language model (top). LLM | Size | Top10 | Bottom10 | Diff.%
---|---|---|---|---
|  | Err.% | Err. % |
Salary |  |
mistral-instruct | 7b | 23.4 | 25.2 | 7.6
llama-2-chat | 7b | 17.4 | 27.2 | 56.2
llama-2-chat | 70b | 14.2 | 20.1 | 41.9
gpt-3.5-turbo | 20b | 15.0 | 23.3 | 55.0
gpt-4 | ? | 13.2 | 23.4 | 77.0
Employer Presence |  |
mistral-instruct | 7b | 240.7 | 330.7 | 37.4
llama-2-chat | 7b | 46.4 | 467.6 | 906.9
llama-2-chat | 70b | 169.5 | 631.9 | 272.8
gpt-3.5-turbo | 20b | 100.1 | 323.0 | 222.9
gpt-4 | ? | 68.6 | 244.2 | 256.0
Commute Duration |  |
mistral-instruct | 7b | 215.0 | 431.7 | 100.8
llama-2-chat | 7b | 155.7 | 268.9 | 72.7
llama-2-chat | 70b | 169.5 | 631.9 | 272.8
gpt-3.5-turbo | 20b | 83.4 | 188.5 | 125.9
gpt-4 | ? | 56.5 | 232.8 | 312.4
Table 2: Performance of language models on the 10 largest and 10 smallest
metropolitan areas. The “Diff.%” column indicates the normalized percentage
increase in error when comparing the bottom 10 to the top 10. i.e. a diff/% of
10 means that the error is 10% higher for the bottom 10 as compared to the top
10.

Experimental results are displayed in Table 1. Nearly all experiments
demonstrate negative and significant Pearson coefficients, indicating that
LLMs tend to achieve better performance in tasks related to larger cities.
Commute duration tasks generally show less stable outcomes, with lower
correlations and a broader range of median errors, whereas salary prediction
tasks consistently present higher correlations and smaller median errors.
Median errors are highest for employer presence predictions, which possibly
relates to its distribution being the most skewed (Figure 1).

##  4 Analysis

###  4.1 Visualizations

Rendering scatterplots of our predictions (Figure 3), we observe the same
trend as in Table 1: all correlations are negative to some extent. In these
visualizations, we observe that the LLaMA series of language models have the
most outliers, with the 70b parameter variant having the most trouble with
larger populations (especially for salary). This suggests that larger model
size does not necessarily lead to more consistent predictions.

###  4.2 Top 10 vs Bottom 10

In Table 2, we compare the performance of the language models on the largest
and smallest 10 metropolitan areas. It is evident that performance is better
for tasks associated with larger regions. While target variable error
magnitudes vary greatly, even the approx. 50% average difference among salary
tasks is quite noticeable. Performance disparities among commute duration and
employer presence are much larger, with up to 9x worse employer presence
performance in the bottom 10 metropolitan areas by llama-2-chat-7b.

##  5 Related Work

Research using the open-source PopQA dataset shows that LLMs exhibit worse
performance at tasks associated with "popular" entities. Mallen et al. (2023)
Additionally, a separate study reveals a correlation between a nation’s GDP-
per-capita and performance on related LLM tasks. Kaplunovich (2023) While both
studies emphasize disparities in LLM task performance based on popularity
and/or geography, neither assesses HR-related task performance across
metropolitan areas.

##  6 Conclusion

Large language models, while powerful, show suboptimal performance in
predicting salaries, commute duration, and employer presence in specific
regions, with this trend worsening in smaller areas. While LLMs seem
unsuitable at generating such job matching data, practitioners should remain
vigilant in mitigating geographic bias as these models undergo further
development and improvement.

## Limitations

Our study’s primary limitation lies in the inherent nature of language models
as probabilistic systems, which leads to inconsistent outputs. For example,
generating llama-2-chat-7b employment presence predictions only yielded 259
valid, non-outlier data points, less than other models. Additionally, our
research was predominantly U.S.-centric, limiting its applicability to other
geographical contexts. Another limitation is our focus on absolute error
metrics, which fail to indicate whether the model systematically overestimates
or underestimates certain variables, like salaries. Addressing these issues in
future research could improve accuracy and applicability.

## References

  * Brown et al. (2020) Tom Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared D Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, et al. 2020.  Language models are few-shot learners.  _Advances in neural information processing systems_ , 33:1877–1901.
  * Chiarain et al. (2023) Neasa-Ni Chiarain, Oisín Nolan, Neimhin-Robinson Gunning, and Madeleine Comtois. 2023.  Filling the slate: examining the contribution llms can make in irish story content generation.  In _Proc. 9th Workshop on Speech and Language Technology in Education (SLaTE)_ , pages 176–181.
  * Jiang et al. (2023) Albert Q Jiang, Alexandre Sablayrolles, Arthur Mensch, Chris Bamford, Devendra Singh Chaplot, Diego de las Casas, Florian Bressand, Gianna Lengyel, Guillaume Lample, Lucile Saulnier, et al. 2023.  Mistral 7b.  _arXiv preprint arXiv:2310.06825_.
  * Kaplunovich (2023) A. Kaplunovich. 2023.  Wealth of nations, wealth of data: How gdp shapes diverse large language models like chatgpt : Interviewing assorted open source generative ai models.  In _2023 IEEE International Conference on Big Data (BigData)_ , pages 4654–4663, Los Alamitos, CA, USA. IEEE Computer Society.
  * Koto et al. (2023) Fajri Koto, Nurul Aisyah, Haonan Li, and Timothy Baldwin. 2023.  Large language models only pass primary school exams in Indonesia: A comprehensive test on IndoMMLU.  In _Proceedings of the 2023 Conference on Empirical Methods in Natural Language Processing_ , pages 12359–12374, Singapore. Association for Computational Linguistics.
  * Mallen et al. (2023) Alex Mallen, Akari Asai, Victor Zhong, Rajarshi Das, Daniel Khashabi, and Hannaneh Hajishirzi. 2023.  When not to trust language models: Investigating effectiveness of parametric and non-parametric memories.  In _Proceedings of the 61st Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 9802–9822, Toronto, Canada. Association for Computational Linguistics.
  * OpenAI et al. (2023) OpenAI, J. Achiam, S. Adler, S. Agarwal, L. Ahmad, I. Akkaya, F. L. Aleman, D. Almeida, J. Altenschmidt, S. Altman, S. Anadkat, R. Avila, I. Babuschkin, S. Balaji, V. Balcom, P. Baltescu, H. Bao, M. Bavarian, and B. Belgum, J. … Zoph. 2023.  Gpt-4 technical report.
  * Pearson (1901) Karl Pearson. 1901.  LIII. on lines and planes of closest fit to systems of points in space.  _The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science_ , 2(11):559–572.
  * Touvron et al. (2023a) Hugo Touvron, Thibaut Lavril, Gautier Izacard, Xavier Martinet, Marie-Anne Lachaux, Timothée Lacroix, Baptiste Rozière, Naman Goyal, Eric Hambro, Faisal Azhar, et al. 2023a.  Llama: Open and efficient foundation language models.  _arXiv preprint arXiv:2302.13971_.
  * Touvron et al. (2023b) Hugo Touvron, Louis Martin, Kevin Stone, Peter Albert, Amjad Almahairi, Yasmine Babaei, Nikolay Bashlykov, Soumya Batra, Prajjwal Bhargava, Shruti Bhosale, et al. 2023b.  Llama 2: Open foundation and fine-tuned chat models.  _arXiv preprint arXiv:2307.09288_.


# Recent advances in the SISSO method and their implementation in the SISSO++
code

Thomas A. R. Purcell<EMAIL_ADDRESS>The NOMAD Laboratory at the
FHI of the Max-Planck-Gesellschaft and IRIS-Adlershof of the Humboldt-
Universität zu Berlin, Faradayweg 4–6, D-14195 Berlin, Germany  Matthias
Scheffler  The NOMAD Laboratory at the FHI of the Max-Planck-Gesellschaft and
IRIS-Adlershof of the Humboldt-Universität zu Berlin, Faradayweg 4–6, D-14195
Berlin, Germany  Luca M. Ghiringhelli<EMAIL_ADDRESS>Physics
Department and IRIS-Adlershof, Humboldt Universität zu Berlin, Zum Großen
Windkanal 2 D-12489 Berlin, Germany  The NOMAD Laboratory at the FHI of the
Max-Planck-Gesellschaft and IRIS-Adlershof of the Humboldt-Universität zu
Berlin, Faradayweg 4–6, D-14195 Berlin, Germany

###### Abstract

Accurate and explainable artificial-intelligence (AI) models are promising
tools for the acceleration of the discovery of new materials, ore new
applications for existing materials. Recently, symbolic regression has become
an increasingly popular tool for explainable AI because it yields models that
are relatively simple analytical descriptions of target properties. Due to its
deterministic nature, the sure-independence screening and sparsifying operator
(SISSO) method is a particularly promising approach for this application. Here
we describe the new advancements of the SISSO algorithm, as implemented into
SISSO++, a C++ code with Python bindings. We introduce a new representation of
the mathematical expressions found by SISSO. This is a first step towards
introducing “grammar” rules into the feature creation step. Importantly, by
introducing a controlled non-linear optimization to the feature creation step
we expand the range of possible descriptors found by the methodology. Finally,
we introduce refinements to the solver algorithms for both regression and
classification, that drastically increase the reliability and efficiency of
SISSO. For all of these improvements to the basic SISSO algorithm, we not only
illustrate their potential impact, but also fully detail how they operate both
mathematically and computationally.

##  I Introduction

Data-centric and artificial intelligence (AI) approaches are becoming a vital
tool for describing physical and chemical properties and processes. The key
advantage of AI is its ability to find correlations between different sets of
properties without the need to know which ones are important before the
analysis. Because of this, AI has become increasingly popular for materials
discovery applications with uses in areas such as thermal transport properties
[1, 2], catalysis [3, 4, 5], phase stability [6], and quantum materials [7].
Despite the success of these methodologies, creating explainable and
physically relevant AI models remains an open challenge in the field [8, 9,
10].

One prevalent set of methods for explainable AI is symbolic regression [11,
12, 13, 14]. Symbolic regression algorithms identify the optimal non-linear,
analytic expressions for a given target property from a set of input features,
i.e., the primary features, that are possibly related to the target [15].
Originally, (stochastic) genetic-programming-based approaches were and still
are used to find these expressions [16, 17, 18, 15], but recently a more
diverse set of solvers have been developed [19, 20, 21, 22, 23, 24]. The sure-
independence screening and sparsifying operator (SISSO) approach combines
symbolic regression with compressed sensing [25, 26, 27, 28], to provide a
deterministic way of finding these analytic expressions. This approach has
been used to describe numerous properties including phase stability [6, 26,
29], catalysis [30], and glass transition temperatures [31]. It has also been
used in a multi-task [25] and hierarchical fashion [28].

The SISSO approach starts with a collection of primary features and
mathematical unary and binary operators (e.g., addition, multiplication, $n$th
root, logarithms, etc.). The first step is the feature-creation step, where a
pool of generated features, is built by exhaustively applying the set of
mathematical operators to the primary features. The algorithm is iteratively
repeated by applying the set of operators to the primary features or the
generated features created at the previous step. The number of iterations in
this feature-creation step is called the rung. The subsequent step is
descriptor identification, i.e., compressed sensing is used to identify the
best $n$-dimensional linear model by performing an $\ell_{0}$-regularized
optimization on a subspace $\mathcal{S}$ of all generated expressions.
$\mathcal{S}$ is selected using sure-independence screening [32], with a
suitable projection score, depending on whether one is solving a regression or
classification problem (see below). For a regression problem this will rank
all generated features according to their Pearson correlation values and
select only the most correlated features, essentially performing a one-
dimensional $\ell_{0}$-regularized linear regression on all features. The
result of the SISSO analysis is a $n-$dimensional descriptor, which is a
vector with components from $\mathcal{S}$. For a regression problem, the SISSO
model is the scalar product of the identified descriptor with the vector of
linear coefficients resulting from the $\ell_{0}$-regularized linear
regression. For a classification problem, the model is given as a set of
hyperplanes that divide the points into classes, that are described by the
scalar product of the identified descriptor, with a set of coefficients found
by linear support vector machines (SVM).

Here, we introduce the new concepts implemented in the recently released
SISSO++ code [27] and detail their implementation. Beyond creating a modular
interface to run SISSO, SISSO++ also improves upon the algorithms in several
aspects, for both the feature creation and descriptor identification steps.
The most important advancement of the code is expressing the features as
binary expression trees, instead of strings, allowing us to recursively define
all aspects of the generated expressions from the primary features. With this
implementation choice, we are able to create a complete description of the
units for each generated feature, as well as an initial representation of its
domain and range. This allows for the creation of grammatically correct
expressions, in terms of consistency of the physical units, and the control of
numerical issues generated by features going out of their physically
meaningful range. In terms of the feature-creation step, we also discuss the
implementation of Parameteric SISSO, which introduces the flexibility of non-
linear parameters together with the operators that are optimized against a
loss function based on the compressed-sensing-based feature selection metrics.
This procedure was used to describe the thermal conductivity of a material in
a recent publication [33]. For the descriptor-identification step, we cover
two components: an improved classification algorithm and the multi-residual
approach. For classification problems, we generalized the algorithm to work
for any problem to an arbitrary dimension, and explicitly include a model
identification via linear SVM. The multi-residual approach, which was
previously used in Ref. [28], introduces further flexibility for the
identification of models with more than one dimension. Here, we provide a in-
depth discussion of its machinery.

##  II Feature Creation

###  II.1 Binary-Expression-Tree Representation of Features

Figure 1: A demonstration of the new representation of the features in the
SISSO++ code. The feature is stored as the root of the tree (represented by
the thick border), the primary features are the leaves, and the rung
corresponds to the height of the tree, i.e. the longest path between each leaf
and the root. The unit, range, expressions, and values of the features are
necessarily stored only for the primary features, with them defined
recursively for all generated expressions.

The biggest advancement of the implementation of SISSO in SISSO++, compared to
the original implementation in Ref. [26], is its modified representation of
the features as binary expression trees, instead of strings. This
representation is illustrated in Figure 1 and easily allows for all aspects of
the generated features to be recursively calculated on the fly from the data
stored in the primary features. For certain applications, it is also possible
to store the data of higher-rung features, to reduce the overall computational
cost of the calculations. The individual features are addressed by the root
node of the binary expression tree, and stored in the code as a shared_ptr
from the C++ standard library. This representation reduces the overall memory
footprint of each calculation, as the individual features only need to be
created once and only copies of shared pointers need to be stored for each new
expression. The remainder of this section will be used to describe the various
aspects of the new representation including a description of the units and
range of the features, as well as how it is used to generate the feature
space.

####  II.1.1 Units

An important upgrade in SISSO++ is its generalized and exact treatment of
units for the expressions. In physics, dimensional analysis is an important
tool when generating physically meaningful expressions, and it is necessary to
include it when using symbolic regression for scientific problems. We
introduced this into SISSO++ by determining the units for each new expression
from the primary features, and explicitly checking to ensure that a new
expression is physically possible. Within the code the Units are implemented
as dictionaries with the key representing the base unit and the value
representing the exponent for each key, e.g., $\mathrm{m/s^{2}}$ would be
stored as $\mathrm{\\{m:1,s:-2\\}}$. Functions exist to transform the Units to
and from strings to more easily represent the information. We then implemented
a multiplication, division, and power functions for these specialized
dictionaries, allowing for the units of the generated features to be derived
recursively following the rules in Table 1. An important caveat is that the
current implementation can not convert between two units for the same physical
quantity units, e.g. between nanometers, picometers, and Bohr radii.

Table 1: How the units are calculated for each operation Operation | Resulting Unit
---|---
$A+B$ | Unit(A)
$A-B$ | Unit(A)
$A\left(B\right)$ | Unit(A) * Unit(B)
$A/B$ | Unit(A) / Unit(B)
$\left|A-B\right|$ | Unit(A)
$\left|A\right|$ | Unit(A)
$\sin\left(A\right)$ | Unitless
$\cos\left(A\right)$ | Unitless
$\exp\left(A\right)$ | Unitless
$\exp\left(-A\right)$ | Unitless
$\log\left(A\right)$ | Unitless
$\left(A\right)^{-1}$ | Unit(A)-1
$\left(A\right)^{2}$ | Unit(A)2
$\left(A\right)^{3}$ | Unit(A)3
$\left(A\right)^{6}$ | Unit(A)6
$\sqrt{A}$ | Unit(A)1/2
$\sqrt[3]{A}$ | Unit(A)1/3

Using this implementation of units, a minimal treatment of dimensional
analysis can be performed in the code. The dimensional analysis focuses on
whether the units are consistent within each expression and for the final
model. For the final model, this check is used to determine the units of the
fitted constants in the linear models at the end, which can take arbitrary
units, and therefore can do any unit conversion natively. The restrictions,
used to reject expressions by dimensional analysis, are outlined in Table 2,
and can be summarized as only allowing addition and subtraction to act on
features of the same units, and all transcendental operations must act on a
unitless quantity. With these two restrictions in place, only physically
possible features can be found, and the choice of units should no longer
affect which features are selected. If one wants to allow for non-physical
expressions to be found, this can be achieved by providing all input primary
features with no units associated to them.

Table 2: Restrictions for each unit, if an operation is not listed there are no restrictions. Operation | Unit Restriction
---|---
$A+B$ | Unit(A) == Unit(B)
$A-B$ | Unit(A) == Unit(B)
$\left|A-B\right|$ | Unit(A) == Unit(B)
$\sin\left(A\right)$ | Unit(A) == $\emptyset$
$\cos\left(A\right)$ | Unit(A) == $\emptyset$
$\exp\left(A\right)$ | Unit(A) == $\emptyset$
$\exp\left(-A\right)$ | Unit(A) == $\emptyset$
$\log\left(A\right)$ | Unit(A) == $\emptyset$

####  II.1.2 Range

Another important advancement of the feature-creation step is the introduction
of ranges for the primary features, which act as a domain for future
operations during feature creation. One of the challenges associated with
symbolic regression, especially with smaller datasets, is that the selected
expressions can sometimes contain discontinuities that are outside of the
training data, but still within the relevant input space for a given problem.
For example, this can lead to an expression taking the logarithm of a negative
number, resulting in an undefined prediction. SISSO++ solves this problem by
including an option for describing the range of a primary feature in standard
mathematical notation, e.g. $\left[\left.0,\infty\right)\right.$, and then
using that to calculate the range for all generated formulas using that
primary feature, following the rule specified in Table 3. In the code, the
ranges are referenced as the Domain because the range for a feature of rung
$n-1$ is the domain for a possible expression of rung $n$ that is using that
feature. While all ranges in Table 3 assume inclusive endpoints, the
implementation can handle both exclusive endpoints and a list of values
explicitly excluded from the range, e.g., point discontinuities inside the
primary features themselves.

Table 3: How the range for each operation is calculated Operation | Resulting Range
---|---
$A+B$ | $\left[\min\left(A\right)+\min\left(B\right),\max\left(A\right)+\max\left(B\right)\right]$
$A-B$ | $\left[\min\left(A\right)-\max\left(B\right),\max\left(A\right)-\min\left(B\right)\right]$
$A\left(B\right)$ | $\left[\min\left(\min\left(A\right)*\min\left(B\right),\min\left(A\right)*\max\left(B\right),\max\left(A\right)*\min\left(B\right),\max\left(A\right)*\max\left(B\right)\right),\right.$
|
$\left.\max\left(\min\left(A\right)*\min\left(B\right),\min\left(A\right)*\max\left(B\right),\max\left(A\right)*\min\left(B\right),\max\left(A\right)*\max\left(B\right)\right)\right]$
$A/B$ | $\left[\mathrm{Range}\left(A\right)*\mathrm{Range}\left(B^{-1}\right)\right]$
$\sin\left(A\right)$ | $\left[-1,1\right]$
$\cos\left(A\right)$ | $\left[-1,1\right]$
$\exp\left(A\right)$ | $\left[\exp\left(\min\left(A\right)\right),\exp\left(\max\left(A\right)\right)\right]$
$\exp\left(-A\right)$ | $\left[\exp\left(-\max\left(A\right)\right),\exp\left(-\min\left(A\right)\right)\right]$
$\log\left(A\right)$ | $\left[\log\left(\min\left(A\right)\right),\log\left(\max\left(A\right)\right)\right]$
$\left(A\right)^{-1}$ | if($0\in\mathrm{Range}\left(A\right)$):
|  $\left(-\infty,0\right)\cup\left(0,\infty\right)$
| else:
|
$\left[\left(\max\left(A\right)\right)^{-1},\left(\min\left(A\right)\right)^{-1}\right]$
$\left(A\right)^{3}$ | $\left[\left(\min\left(A\right)\right)^{3},\left(\max\left(A\right)\right)^{3}\right]$
$\sqrt{A}$ | $\left[\sqrt{\min\left(A\right)},\sqrt{\max\left(A\right)}\right]$
$\sqrt[3]{A}$ | $\left[\sqrt[3]{\min\left(A\right)},\sqrt[3]{\max\left(A\right)}\right]$
$\left|A\right|$ | $\left[\max\left(0,\min\left(A\right)\right),\max\left(\left|\max\left(A\right)\right|,\left|\min\left(A\right)\right|\right)\right]$
$\left(A\right)^{2}$ | $\left[\max\left(0,\min\left(A\right)\right)^{2},\max\left(\left|\max\left(A\right)\right|,\left|\min\left(A\right)\right|\right)^{2}\right]$
$\left(A\right)^{6}$ | $\left[\max\left(0,\min\left(A\right)\right)^{6},\max\left(\left|\max\left(A\right)\right|,\left|\min\left(A\right)\right|\right)^{6}\right]$
$\left|A-B\right|$ | $\left[\max\left(0,\min\left(A\right)-\max\left(B\right)\right),\right.$
|
$\left.\max\left(\left|\max\left(\max\left(A\right)-\min\left(B\right)\right)\right|,\left|\min\left(\max\left(A\right)-\min\left(B\right)\right)\right|\right)\right]$

Table 4 lists the cases where the range of a feature is used to prevent a new
expression from being generated. In all cases, this prevents an operation from
occurring where a mathematical operation would be not defined, such as taking
the square-root of a negative number. In cases where the range of values for a
primary feature is not defined, then these checks are not performed and the
original assumption that all operations are safe is used.

Table 4: Domain restrictions for each operation, if an operation is not listed there are no restrictions. Operation | Domain Restriction
---|---
$A/B$ | $0\notin\mathrm{Range\left(B\right)}$
$\left(A\right)^{-1}$ | $0\notin\mathrm{Range\left(A\right)}$
$\log\left(A\right)$ | $\min\left(\mathrm{Range}\left(A\right)\right)>0$
$\sqrt{A}$ | $\min\left(\mathrm{Range}\left(A\right)\right)\geq 0$

###  II.2 Parametric SISSO

Figure 2: A graphical representation of the effect of the parameterization
depth. If $P_{d}=1$ then only the $\mathrm{sin}$ operator gets parameterized,
while if $P_{d}=2$ both operations are parameterized

Parametric SISSO extends the feature creation step of SISSO to automatically
include scale and bias terms for each operation, as used in Purcell et
al.[33]. For a general operator,
$\hat{h}\left(\mathbf{x}\right)\in\hat{\mathcal{H}}$, with a set of scale and
bias parameters, $\hat{\mathcal{P}}$, the parameterization scheme updates the
operator to be,

$\hat{h}\left(\mathbf{x}\right)\to\hat{h}^{\hat{\mathcal{P}}}\left(\alpha_{1}\mathbf{x}+\beta_{1}\right),$
(1)

where $\alpha_{1}$ is the scale parameter, $\beta_{1}$ is the bias term, and
$\mathbf{x}$ is a vector containing all input data. These new operators can
then be used to create a new feature,
$\hat{\phi}^{\hat{\mathcal{P}}}\left(\mathbf{x}\right)$, as is normally done
in SISSO, where each feature has its own set of parameters. However, this does
introduce a new hyperparameter for the feature creation step, the
parameterization depth, $P_{d}$, which defines the level at which the
parameterization occurs in the binary expression tree. This is best described
in Figure 2, where $P_{d}$ controls which operations are included in the
parameterization. In this example, if $\beta_{2}$ was previously set by
another optimization and $P_{d}=2$, then that previous value will be ignored
for this feature only; however, if $P_{d}=1$, then the existing $\beta_{2}$
will be preserved. In order to avoid linear dependencies between different
operations and the constants set during linear regression, some of the scale
and bias terms are set to one and zero, respectively, for a summary of this
for all operators see Table 5. It is important to note that for the log
operator, $\alpha$ is always set to $\pm 1$ to avoid linear dependencies with
other parameters. Although this does leave a unit dependency, it can be
removed with

$\ln\left(x+\beta\right)\rightarrow\ln\left(\alpha_{\mathrm{unit}}\left(x+\beta\right)\right)-\ln\left(\alpha_{\mathrm{unit}}\right),$
(2)

where $\alpha_{\mathrm{unit}}$ is the unit conversion factor.

Figure 3: A comparison of the expressions found non-parametric (a, c, and e) and parametric SISSO (b, d, and f) for a Lorentzian (a, b, e, and f) and sin (c, d) function. Blue dots represent the training data, and the red line represents the expressions found by SISSO. The parameterization scheme either finds the correct (b and d) or better model (f) than the non-parametric functions, even when the high noise or bad initial guess of the parameters leads to a non-optimized solution. Table 5: The free parameters and updated equations for each operation used in the updated feature creation step of SISSO Operation | Parameterized Operation
---|---
$A+B$ | $A+\alpha B$
$A-B$ | $A-\alpha B$
$A\left(B\right)$ | $A\left(B+\beta\right)$
$A/B$ | $A/\left(B+\beta\right)$
$\left|A-B\right|$ | $\left|A-\left(\alpha B+\beta\right)\right|$
$\left|A\right|$ | $\left|A+\beta\right|$
$\sin\left(A\right)$ | $\sin\left(\alpha A+\beta\right)$
$\cos\left(A\right)$ | $\cos\left(\alpha A+\beta\right)$
$\exp\left(A\right)$ |  $\exp\left(\alpha A\right)$111$\alpha>0$
$\exp\left(-A\right)$ |  $\exp\left(-\alpha A\right)$11footnotemark: 1
$\log\left(A\right)$ |  $\log\left(A+\beta\right)$222$\alpha$ can only be 1.0 or -1.0
$\left(A\right)^{-1}$ | $\left(A+\beta\right)^{-1}$
$\left(A\right)^{2}$ | $\left(A+\beta\right)^{2}$
$\left(A\right)^{3}$ | $\left(A+\beta\right)^{3}$
$\left(A\right)^{6}$ | $\left(A+\beta\right)^{6}$
$\sqrt{A}$ |  $\sqrt{\alpha A+\beta}$22footnotemark: 2
$\sqrt[3]{A}$ | $\sqrt[3]{A+\beta}$

Once $\hat{\phi}^{\hat{\mathcal{P}}}\left(\mathbf{x}\right)$ is defined, all
parameters $\hat{p}\in\hat{\mathcal{P}}$ are optimized using the non-linear
optimization library NLopt [34]. We use the Cauchy loss function as the
objective for the optimization

$\displaystyle\underset{\hat{\mathcal{P}}}{\mathrm{min}}$
$\displaystyle\,f\left(\mathbf{P},\hat{\phi}^{\hat{\mathcal{P}}}\right)$ (3a)
$\displaystyle f\left(\mathbf{P},\hat{\phi}^{\hat{\mathcal{P}}}\right)$
$\displaystyle=\sum_{i}^{n_{\mathrm{samp}}}\frac{c^{2}}{n_{\mathrm{samp}}}\log\left(1+\left(\frac{P_{i}-\hat{\phi}^{\hat{\mathcal{P}}}\left(\mathbf{x}_{i}\right)}{c}\right)^{2}\right),$
(3b)

where $\mathbf{P}$ is a property vector, $c$ is a scaling factor set to 0.5
for all calculations and $n_{\mathrm{samp}}$ is the number of samples. We
chose to use the Cauchy loss function over the mean square error to make the
non-linear optimization more robust against outliers in the dataset. Because
Equation 3b is not scale or bias invariant, additional external parameters
$\alpha_{ext}$ and $\beta_{ext}$ are introduced to respectively account for
these effects. For the case of multi-task SISSO [25], each task has its own
external bias and scale parameters to account for the individual linear
regression solutions. As an example for the feature illustrated in Figure 2
($P_{d}=2$), the function that is optimized would be

$\hat{\phi}^{\hat{\mathcal{P}}}\left(\mathbf{x}\right)=\alpha_{ext}\sin\left(\alpha_{1}\frac{x_{1}}{x_{2}+\beta_{2}}+\beta_{1}\right)+\beta_{ext}.$
(4)

To initialize the parameters in $\hat{\mathcal{P}}$, we set all internal
$\alpha$ and $\beta$ terms are set to 1.0 and 0.0, respectively, and
$\alpha_{ext}$ and $\beta_{ext}$ are set to the solution of the least squares
regression problem for each task. In some cases, $\beta$ can be set to a non-
zero value if leaving it at zero would include values outside the Domain of
the operator. In these cases, $\beta$ is set to
$\min\left(\mathbf{x}\right)+10^{-10}$.

Figure 4: Illustration of how the feature space of SISSO is created. In this
example the user selects two primary features $\omega$ (purple) and $t$
(orange) and three operators sin (blue), multiplication (green), and division
(red). SISSO then builds up a more complicated expression space by applying
the operations onto the existing features by increasing the height of the
binary expression trees. Throughout this process the units and ranges of each
of the operations are respected.

Each optimization follows a two or three step process outlined here. First a
local optimization is performed to find the local minimum associated with the
initial parameters. Once at a local minimum an optional global optimization is
performed to find any minima that are better than the one initially found. For
these first two steps the parameters are optimized to a relative tolerance of
$10^{-3}$ and $10^{-2}$ respectively, with a maximum of five thousand function
evaluations. Finally, a more accurate local optimization is done to a relative
tolerance of $10^{-6}$ to find the best parameter set. For this final
optimization, ten thousand function evaluations are allowed. Additionally, for
both the initial and global optimization steps, the parameters are bounded to
be in a range between -100 and 100 to improve the efficiency of the
optimization, but this restriction is removed for the final optimization. For
all local optimizations, the subplex algorithm [35], a faster and more robust
variant of the Nelder-Mead Simplex method [36], is used. The Improved
Stochastic Ranking Evolution Strategy algorithm [37] is used for all global
optimizations. Once optimized, only the internal $\alpha$ and $\beta$
parameters are stored in $\hat{\mathcal{P}}$.

Figure 3 illustrates the power of the new parameterization scheme. For both
toy problems represented by analytic Lorentzian and sin functions with some
white noise, the non-parametric version of SISSO can not find accurate models
for the equations as it can not address the non-linearities properly. By using
this new parameterization scheme, SISSO is now able to accurately find the
models as shown in Figure 3c and d. However, it is important to note that the
more powerful featurization comes at the cost of a significantly increased
time to generate the feature space, as the parameterization becomes the
bottleneck for the calculations. Additionally, there can be cases where the
parameterization scheme does not find an optimal solution because of too much
noise or the optimal parameters being too far away from the initial guesses,
as shown in Figure 3e and f.

###  II.3 Building the Complete Feature Set

With all of the new aspects of the feature representation in place, SISSO++
has a fully parallelized feature set construction that uses a combination of
threads and MPI ranks, for efficient feature set construction. The basic
process of creating new features is illustrated in Figure 4, where each new
rung builds on top of existing features by adding a new operation on top of
existing binary expression trees for the previous rung. Throughout this
process all checks are done to ensure that the units are correct and the
domains for each new operation are respected. Additionally, the code checks
for invalid values, e.g., NaN or inf, and some basic simplifications for all
features, e.g., features like $\frac{t\omega}{\omega}$ are rejected. The
operators are separated into parameterized and non-parameterized versions of
each other to allow for the optional use of the parametric SISSO concepts.

##  III Descriptor Identification

###  III.1 Linear Programming Implementation for Classification Problems

One of the largest updates to the SISSO methodology is the new, generalized
approach for solving classification problems. In previous implementations,
when finding a classification scheme, SISSO would explicitly build the convex
hull, and then calculate the number of points inside the overlap region
between different classes and either the normalized overlap volume or
separation distance to find the optimal solution. While this works for two
dimensions, finding the overlap volume or separation distance becomes
intractable for three or more dimensions, and even defining the convex hull
becomes intractable for four or more dimensional classification. SISSO++
replaces these conditions with an algorithm that determines the number of
points inside the convex hull overlap region using linear programming, and
explicitly creates a model using linear SVM. The linear programming algorithm
checks for the feasibility of

$\displaystyle\min 0$ (5) $\displaystyle\mathrm{s.t.}\sum_{i\in
I}^{\min}\alpha_{i}x_{i}=x_{j},\quad\sum_{i\in
I}\alpha_{i}=1,\quad\alpha_{i}\geq 0,\ \forall i\in I\mathrm{,}$ (6)

where $x_{i}$ is the $i^{th}$ point inside the set of all points of a class
$I$, $\alpha_{i}$ is the coefficient for $x_{i}$, and $x_{j}$ is the point to
check if it is inside the convex hull. The above problem is only feasible if
and only if $x_{j}$ lies inside the set of points, $I$, representing a class
in the problem. Here we are optimizing a zero function because the actual
solution to this optimization does not matter, only that such a solution can
be found, i.e. the constraints can be fulfilled, is important. The feasibility
and linear programming problem is defined using the Coin-CLP library [38].
Once the number of points in the overlap region is determined, a linear SVM
model is calculated for the best candidates, and used as the new tie-breaking
procedure. The first tiebreaker is the number of misclassified points by the
SVM model and the second one is the margin distance. The SVM model is
calculated using libsvm [39].

Figure 5: A demonstration of the classification algorithm. A set of 1000
points of three features, $x_{0}$, $x_{1}$, and $x_{2}$, are sampled from a
Gaussian distribution with a standard deviation of 0.5 and centered at the
origin. The set is separated into 2 classes (the red and light blue circles)
by the plane $x_{0}+x_{1}+x_{2}=0$ (green surface). a) All points where
$\left|x_{1}+x_{2}+x_{3}\right|<0.4$ are moved to a random point from that has
a distance of $0.2\sqrt{2}$ to $0.35\sqrt{2}$ from the plane within each of
their own classes, and b) the original data stored as $x_{3}$, $x_{4}$, and
$x_{5}$. The updated classification algorithm correctly determines that
$x_{0}$, $x_{1}$, and $x_{2}$ is the superior classifier, while for the
original definitions of only the convex overlap region for three and more
dimensions they would be considered equally good. The projections are shown
with an elevation angle of 7.5 degrees and an azimuthal angle of 145 deg.

Figure 5 demonstrates the capabilities of the new classification algorithm. In
this problem, we randomly sample a Gaussian distribution with a standard
deviation of 0.5 and centered at the origin to generate one thousand samples
for three features, $x_{0}$, $x_{1}$, and $x_{2}$. These three features and
the copied and relabeled to $x_{3}$, $x_{4}$, and $x_{5}$ The samples are then
separated into two classes based on if $\left(x_{1}+x_{2}+x_{3}\right)$ is
greater than (light blue) or less than (red) zero. For $x_{0}$, $x_{1}$, and
$x_{2}$ we move all points where $\left|x_{1}+x_{2}+x_{3}\right|<0.4$ to a
random point further away from the dividing plane within each of their classes
to ensure a noticeable margin of separation. For simplicity, we use only the
six primary features shown on the axes of Figure 5a and b with no generated
expressions; however, the rung 2 feature of $\left(x_{0}+x_{1}\right)+x_{2}$
would be able to completely separate the classes in one dimension. Using the
updated algorithm SISSO can now easily identify that the set
$\left\\{x_{0},x_{1},x_{2}\right\\}$ is the better classifier than the set
$\left\\{x_{3},x_{4},x_{5}\right\\}$ as the tiebreakers are defined for all
dimensions and not just the first two. More importantly, SISSO can now provide
the $n$-dimensional dividing plane found by linear-SVM creating an actual
classifier automatically here shown by the green plane. It is important to
note that the same procedure can be done to an arbitrary dimension, but
visualization becomes impractical above three.

###  III.2 Multiple Residuals

Figure 6: An illustration of how tracking multiple residuals can improve the
performance of SISSO. A three dimensional problem space (three training
samples) defined by $\mathbf{e_{1}}$, $\mathbf{e_{2}}$, and $\mathbf{e_{3}}$
with five feature vectors $\mathbf{F_{1}}$ (gray), $\mathbf{F_{2}}$ (blue),
$\mathbf{F_{3}}$ (red), $\mathbf{F_{4}}$ (brown), and $\mathbf{F_{5}}$
(turquoise) for a property vector, $\mathbf{P}$ (purple). For this example the
size of the SIS subspace, $n_{sis}$, is two. The residuals for the best
($\mathbf{\Delta_{1}}$, gray) and second best($\mathbf{\Delta_{2}}$, blue)
one-dimensional models are shown as dashed lines. If the number of residuals,
$n_{res}$, is one, then only $\mathbf{\Delta_{1}}$ is used (d, e) and
$\mathbf{F_{3}}$ will not be selected in the second SIS step. This means the
best two dimensional model is not found. However, if $n_{res}=2$ (b, c), then
$\mathbf{F_{3}}$ is selected and the best two dimensional model, a combination
of $\mathbf{F_{2}}$ and $\mathbf{F_{3}}$, can be found in the second
$\ell_{0}$ step.

The second advancement to the descriptor-identification step of the SISSO
algorithm is the introduction of a multiple-residuals approach to select the
expressions (features) for models with a dimension higher than one. In the
original SISSO algorithm [40], the residual of the previously found model,
$\vec{\Delta}_{D-1}^{0}$, i.e., the difference between the vector storing the
values of the property for each sample, $\vec{P}$, and the estimates predicted
by the ($D\\!-\\!1$)-dimensional model
($\vec{\Delta}_{D-1}^{0}=\vec{P}_{D-1}-\vec{P}$), is used to calculate the
projection score of the candidate features during the SIS step for the best
$D$-dimensional model,
$s_{j}^{0}=R^{2}\left(\vec{\Delta}_{D-1}^{0},\vec{d}_{j}\right)$. Here, $R$ is
the Pearson correlation coefficient, representing a regression problem, and
$j$ corresponds to each expression generated during the feature-creation step
of SISSO. In SISSO++[27], we extend the residual definition and use the best
$r$ residuals to calculate the projection score:
$\mathrm{max}\left(s_{j}^{0},s_{j}^{1},\ldots,s_{j}^{r-1}\right)$. The
multiple-residual concept generalizes the descriptor identification step of
SISSO, by using information from an ensemble of models to determine which
features to add to the selected subspace. The value of $r$ for a calculation
is set as any hyperparmeter via cross-validation

This process is illustrated in Fig. 6 for a three-dimensional problem space
(three training samples) defined by $\mathbf{e_{1}}$, $\mathbf{e_{2}}$, and
$\mathbf{e_{3}}$ with five candidate features
$\mathbf{F_{1}},\cdots,\mathbf{F_{5}}$ to describe the property $\mathbf{P}$.
In principle, there can be an arbitrary number of feature vectors, but for
clarity we only show five. If only a single residual is used, then the
selected two dimensional model will comprise of $\mathbf{F_{1}}$ and
$\mathbf{F_{4}}$, as $\mathbf{F_{3}}$ is never selected because it has the
smallest projection score from the residual of the model found using $F_{1}$.
However, because $\mathbf{F_{2}}$ has a component along $\mathbf{e_{2}}$, a
better two-dimensional model consisting of a linear combination of
$\mathbf{F_{2}}$ and $\mathbf{F_{3}}$ exists, despite $\mathbf{F_{2}}$ not
being the most correlated feature to $\mathbf{P}$. When going to higher
dimensional feature spaces, it becomes more likely that the feature vectors
similarly correlated to the property contain orthogonal information, thus the
need for using multiple residuals in SISSO.

Figure 7: The training RMSE of a two dimensional model against the size of the
SIS subspace size for $y=5.5+0.4158d_{1}-0.0974d_{2}+\delta$, where
$d_{1}=x_{0}^{2}\sqrt[3]{x_{1}}$, $d_{2}=\left|x_{2}^{3}\right|$, and $\delta$
is a Gaussian white noise term pulled from a distribution with a standard
deviation of 0.05. The solid dark blue line is the error when learning using
one residual, the light blue dashed line is the error when learning using
twenty five residuals, and the red solid line is the error when learning with
fifty residuals.

In order to demonstrate the effect of learning over multiple residuals and get
an estimate of the optimal number of residuals and size of the SIS subspace
($n_{\text{sis}}$), we plot the training RMSE for the two dimensional models
for the function$y=5.5+0.4158d_{1}-0.0974d_{2}+\delta$, where
$d_{1}=x_{0}^{2}\sqrt[3]{x_{1}}$, $d_{2}=\left|x_{2}^{3}\right|$, and $\delta$
is a Gaussian white noise term pulled from a distribution with a standard
deviation of 0.05 in Fig. 7. For this problem the best one dimensional
descriptor is $\frac{1}{\sqrt[3]{x_{3}}}$, with $x_{3}$ being explicitly set
to $\left(y+\Delta\right)^{-3}$ and $\Delta$ is Gaussian white noise with a
standard deviation of 20.0. This primary feature was explicitly added in order
to ensure that the best one-dimensional model would not be either $d_{1}$ or
$d_{2}$ in this synthetic problem. Because of this, when using a single
residual the SIS subspace size has to be increased to over 400, before $y$ can
be reproduced by SISSO. However, by increasing the number of residuals to 50,
SISSO can now find which features are most correlated to the residual of
$d_{1}$ and it immediately finds $y$. In a recent paper we published, we
further demonstrate this approaches effectiveness for learning models of the
bulk modulus of cubic perovskites, [28].

##  IV Conclusions

In this paper, we described recently developed improvements to the SISSO
method and their implementation in SISSO++ code, both in terms of their
mathematical and computational details, which constitute a large leap forwards
in terms of the expressivity of the SISSO method. Utilizing these features
provides greater flexibility and control over the expressions found by SISSO,
and acts as a start to introducing “grammatical” rules into SISSO and symbolic
regression. In particular, concepts such as the units and ranges of the
formula could be extended to prune the search space of possible expressions
for the final models. We have also described the implementation of Parameteric
SISSO, which considerably opens up the range of possible expressions found by
SISSO. Finally, we discussed two improvements related to the SISSO solver,
i.e., a linear programming implementation for the classification problems and
the multiple-residuals technique, both providing extended flexibility in the
descriptors and models found by SISSO.

##  V Acknowledgements

T.A.R.P. thanks Christian Carbogno for valuable discussions related to the
parametric SISSO scheme and proof reading those parts of the manuscript. T.P.
thanks Lucas Foppa for discussions related to the multi-residual approach and
proof reading those parts of the manuscript. This work was funded by the NOMAD
Center of Excellence (European Union’s Horizon 2020 research and innovation
program, grant agreement Nº 951786), the ERC Advanced Grant TEC1p (European
Research Council, grant agreement Nº 740233), BigMax (the Max Planck Society’s
Research Network on Big-Data-Driven Materials-Science), and the project
FAIRmat (FAIR Data Infrastructure for Condensed-Matter Physics and the
Chemical Physics of Solids, German Research Foundation, project Nº 460197019).
T.P. would like to thank the Alexander von Humboldt (AvH) Foundation for their
support through the AvH Postdoctoral Fellowship Program.

##  VI Conflict of Interest Statement

The authors have no conflicts to disclose.

##  VII Author Contributions

T.A.R.P. implemented all methods and performed all calculations. T.A.R.P.
ideated all methods with assistance from LMG. MS and LMG supervised the
project. All authors wrote the manuscript.

##  VIII Data Availability Statement

The data that support the findings of this study and all scripts used to
generate the figures are openly available in FigShare at DOI:
10.6084/m9.figshare.22691857.

## References

  * [1] T. Zhu et al.,  Energy Environ. Sci. 14, 3559 (2021).
  * [2] S. A. Miller et al.,  Chem. Mater. 29, 2494 (2017).
  * [3] K. Tran and Z. W. Ulissi,  Nature Catalysis 1, 696 (2018).
  * [4] L. Foppa et al.,  MRS Bulletin 46, 1016 (2021).
  * [5] L. Foppa et al.,  ACS Catalysis 12, 2223 (2022).
  * [6] C. J. Bartel et al.,  Science Advances 5 (2019).
  * [7] V. Stanev, K. Choudhary, A. G. Kusne, J. Paglione, and I. Takeuchi,  Commun. Mater. 2, 105 (2021).
  * [8] P. P. Angelov, E. A. Soares, R. Jiang, N. I. Arnold, and P. M. Atkinson,  WIREs Data Mining and Knowledge Discovery 11, e1424 (2021).
  * [9] D. Gunning et al.,  Science Robotics 4, eaay7120 (2019).
  * [10] A. Das and P. Rad,  arXiv preprint arXiv:2006.11371 (2020).
  * [11] F. Xu et al.,  Explainable AI: A Brief Survey on History, Research Areas, Approaches and Challenges,  in Lect. Notes Comput. Sci. (including Subser. Lect. Notes Artif. Intell. Lect. Notes Bioinformatics), volume 11839 LNAI, pages 563–574, Springer, 2019.
  * [12] G. S. I. Aldeia and F. O. De França,  Measuring feature importance of symbolic regression models using partial effects,  in GECCO 2021 - Proc. 2021 Genet. Evol. Comput. Conf., pages 750–758, 2021.
  * [13] A. Holzinger, A. Saranti, C. Molnar, P. Biecek, and W. Samek,  Explainable AI Methods - A Brief Overview,  in Lect. Notes Comput. Sci. (including Subser. Lect. Notes Artif. Intell. Lect. Notes Bioinformatics), volume 13200 LNAI, pages 13–38, Springer Science and Business Media Deutschland GmbH, 2022.
  * [14] Z. Li, J. Ji, and Y. Zhang,  arXiv: 2111.12210 (2021).
  * [15] Y. Wang, N. Wagner, and J. M. Rondinelli,  MRS Communications 9, 793 (2019).
  * [16] J. R. Koza,  Statistics and Computing 4, 87 (1994).
  * [17] T. Mueller, E. Johlin, and J. C. Grossman,  Physical Review B 89, 115202 (2014).
  * [18] F. Yuan and T. Mueller,  Scientific Reports 7, 17594 (2017).
  * [19] S.-M. Udrescu and M. Tegmark,  arXiv preprint:1905.11481 (2019).
  * [20] S. Kim et al.,  IEEE Transactions on Neural Networks and Learning Systems 32, 4166 (2021).
  * [21] M. D. Cranmer, R. Xu, P. Battaglia, and S. Ho,  Learning symbolic physics with graph networks, 2019.
  * [22] M. Valipour, B. You, M. Panju, and A. Ghodsi,  Symbolicgpt: A generative transformer model for symbolic regression, 2021\.
  * [23] B. K. Petersen et al.,  Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients, 2021.
  * [24] W. Tenachi, R. Ibata, and F. I. Diakogiannis,  Deep symbolic regression for physics guided by units constraints: toward the automated discovery of physical laws, 2023.
  * [25] R. Ouyang, E. Ahmetcik, C. Carbogno, M. Scheffler, and L. M. Ghiringhelli,  J. Phys. Mater. 2, 24002 (2019).
  * [26] R. Ouyang, S. Curtarolo, E. Ahmetcik, M. Scheffler, and L. M. Ghiringhelli,  Phys. Rev. Mater. 2, 83802 (2018).
  * [27] T. A. R. Purcell, M. Scheffler, C. Carbogno, and L. M. Ghiringhelli,  J. Open Source Softw. 7, 3960 (2022).
  * [28] L. Foppa, T. A. Purcell, S. V. Levchenko, M. Scheffler, and L. M. Ghringhelli,  Phys. Rev. Lett. 129, 55301 (2022).
  * [29] G. R. Schleder, C. M. Acosta, and A. Fazzio,  ACS Appl. Mater. Interfaces 12, 20149 (2020).
  * [30] Z.-K. Han et al.,  Nat. Commun. 12, 1833 (2021).
  * [31] G. Pilania, C. N. Iverson, T. Lookman, and B. L. Marrone,  J. Chem. Inf. Model. 59, 5013 (2019).
  * [32] J. Fan and J. Lv,  J. R. Stat. Soc. Ser. B Stat. Methodol. 70, 849 (2008).
  * [33] T. A. R. Purcell, M. Scheffler, L. M. Ghiringhelli, and C. Carbogno,  submitted to npj Comput. Mater. (2022).
  * [34] S. G. Johnson,  The NLopt nonlinear-optimization package,  http://github.com/stevengj/nlopt, 2021.
  * [35] T. H. Rowan,  Functional stability analysis of numerical algorithms,  PhD thesis, University of Texas at Austin, 1990.
  * [36] J. A. Nelder and R. Mead,  Comput. J. 7, 308 (1965).
  * [37] T. Runarsson and X. Yao,  IEEE Trans. Syst. Man Cybern. Part C (Applications Rev. 35, 233 (2005).
  * [38] J. J. Forrest et al.,  coin-or/clp: Version 1.17.6.
  * [39] C.-C. Chang and C.-J. Lin,  ACM Transactions on Intelligent Systems and Technology 2, 27:1 (2011),  Software available at http://www.csie.ntu.edu.tw/~cjlin/libsvm.
  * [40] R. Ouyang, S. Curtarolo, E. Ahmetcik, M. Scheffler, and L. M. Ghiringhelli,  Physical Review Materials 2, 083802 (2018).