SymbolicRegression.jl is a Julia library for discovering symbolic models to solve optimization problems, especially for regression. It is highly optimized and permits multi-node distributed searches. The library is now very modular and extensible, supporting expressions with custom operators, and regressing on custom types. See the docs for some examples: Home · SymbolicRegression.jl.
It’s hard to believe it’s nearly been 4 years since my original post, so I wanted to give an update to mark the 1.0.0 release! As a summary of the most recent changes, I made this video:
Here’s a bit of a quickstart on using SymbolicRegression. I’ll show the release notes at the end.
The easiest way to use SymbolicRegression.jl is with MLJ.jl. Here’s an example:
import SymbolicRegression: SRRegressor
import MLJ: machine, fit!, predict, report
# Dataset with two named features:
X = (a = rand(500), b = rand(500))
# and one target:
y = @. 2 * cos(X.a * 23.5) - X.b ^ 2
# with some noise:
y = y .+ randn(500) .* 1e-3
model = SRRegressor(
niterations=50,
binary_operators=[+, -, *],
unary_operators=[cos],
)
Now, let’s create and train this model on our data:
mach = machine(model, X, y)
fit!(mach)
Let’s look at the expressions discovered:
report(mach)
Finally, we can make predictions with the expressions on new data:
predict(mach, X)
This will make predictions using the expression selected by model.selection_method, which by default is a mix of accuracy and complexity.
You can override this selection and select an equation from the Pareto front manually with:
predict(mach, (data=X, idx=2))
where here we choose to evaluate the second equation.
For fitting multiple outputs, one can use MultitargetSRRegressor (and pass an array of indices to idx in predict for selecting specific equations). For a full list of options available to each regressor, see the API page.
1.0.0 Release Notes
Now, in the 1.0.0 there are a few major changes compared with the initial release. I give an overview of these here, and then add more detail below.
- Changed the core expression type from
Node{T}→Expression{T,Node{T},Metadata{...}}- This gives us new features, improves user hackability, and greatly improves ergonomics!
- Created “Template Expressions”, for fitting expressions under a user-specified functional form (
TemplateExpression <: AbstractExpression)- Template expressions are quite flexible: they are a meta-expression that wraps multiple other expressions, and combines them using a user-specified function.
- This enables vector expressions - in other words, you can learn multiple components of a vector, simultaneously, with a single expression! Or more generally, you can learn expressions onto any Julia struct.
- (Note that this still does not permit learning using non-scalar operators, though we are working on that!)
- Template expressions also make use of colored strings to represent each part in the printout, to improve readability.
- Created “Parametric Expressions”, for custom functional forms with per-class parameters: (
ParametricExpression <: AbstractExpression)- This lets you fit expressions that act as models of multiple systems, with per-system parameters!
- Introduced a variety of new abstractions for user extensibility (and to support new research on symbolic regression!)
AbstractExpression, for increased flexibility in custom expression types.mutate!andAbstractMutationWeights, for user-defined mutation operators.AbstractSearchState, for holding custom metadata during searches.AbstractOptionsandAbstractRuntimeOptions, for customizing pretty much everything else in the library via multiple dispatch. Please make an issue/PR if you would like any particular internal functions be declaredpublicto enable stability across versions for your tool.- Many of these were motivated to modularize the implementation of LaSR, an LLM-guided version of SymbolicRegression.jl, so it can sit as a modular layer on top of SymbolicRegression.jl.
- Added TensorBoardLogger.jl and other logging integrations via
SRLogger - Support for Zygote.jl and Enzyme.jl within the constant optimizer, specified using the
autodiff_backendoption - Other changes:
- Fundamental improvements to the underlying evolutionary algorithm
- New mutation operators introduced,
swap_operandsandrotate_tree– both of which seem to help kick the evolution out of local optima. - New hyperparameter defaults created, based on a Pareto front volume calculation, rather than simply accuracy of the best expression.
- New mutation operators introduced,
- Changed output file handling
- Major refactoring of the codebase to improve readability and modularity
- Identified and fixed a major internal bug involving unexpected aliasing produced by the crossover operator
- Segmentation faults caused by this are a likely culprit for some crashes reported during multi-day multi-node searches.
- Introduced a new test for aliasing throughout the entire search state to prevent this from happening again.
- Improved progress bar and StyledStrings integration.
- Julia 1.10 is now the minimum supported Julia version.
- Other small features
- Also see the “Update Guide” below for more details on upgrading.
- New URL: https://ai.damtp.cam.ac.uk/symbolicregression
- Fundamental improvements to the underlying evolutionary algorithm
Note that some of these features were recently introduced in patch releases since they were backwards compatible. I am noting them here for visibility.
1. Changed the core expression type from Node{T} → Expression{T,Node{T},...}
This is a breaking change in the format of expressions returned by SymbolicRegression. Now, instead of returning a Node{T}, SymbolicRegression will return a Expression{T,Node{T},...} (both from equation_search and from report(mach).equations). This type is much more convenient and high-level than the Node type, as it includes metadata relevant for the node, such as the operators and variable names.
This means you can reliably do things like:
using SymbolicRegression: Options, Expression, Node
options = Options(binary_operators=[+, -, *, /], unary_operators=[cos, exp, sin])
operators = options.operators
variable_names = ["x1", "x2", "x3"]
x1, x2, x3 = [Expression(Node(Float64; feature=i); operators, variable_names) for i=1:3]
## Use the operators directly!
tree = cos(x1 - 3.2 * x2) - x1 * x1
You can then do operations with this tree, without needing to track operators:
println(tree) # Looks up the right operators based on internal metadata
X = randn(3, 100)
tree(X) # Call directly!
tree'(X) # gradients of expression
Each time you use an operator on or between two Expressions that include the operator in its list, it will look up the right enum index, and create the correct Node, and then return a new Expression.
You can access the tree with get_tree (guaranteed to return a Node), or get_contents – which returns the full info of an AbstractExpression, which might contain multiple expressions (which get stitched together when calling get_tree).
2. Created “Template Expressions”, for fitting expressions under a user-specified functional form (TemplateExpression <: AbstractExpression)
Template Expressions allow users to define symbolic expressions with a fixed structure, combining multiple sub-expressions under user-specified constraints. This is particularly useful for symbolic regression tasks where domain-specific knowledge or constraints must be imposed on the model’s structure.
This also lets you fit vector expressions using SymbolicRegression.jl, where vector components can also be shared!
A TemplateExpression is constructed by specifying:
- A named tuple of sub-expressions (e.g.,
(; f=x1 - x2 * x2, g=1.5 * x3)). - A structure function that defines how these sub-expressions are combined in different contexts.
For example, you can create a TemplateExpression that enforces the constraint: sin(f(x1, x2)) + g(x3)^2 - where we evolve f and g simultaneously.
To do this, we first describe the structure using TemplateStructure that takes a single closure function that maps a named tuple of ComposableExpression expressions and a tuple of features:
using SymbolicRegression
structure = TemplateStructure{(:f, :g)}(
((; f, g), (x1, x2, x3)) -> sin(f(x1, x2)) + g(x3)^2
)
This defines how the TemplateExpression should be evaluated numerically on a given input.
The number of arguments allowed by each expression object is inferred using this closure, though it can also be passed explicitly with the num_features kwarg.
operators = Options(binary_operators=(+, -, *, /)).operators
variable_names = ["x1", "x2", "x3"]
x1 = ComposableExpression(Node{Float64}(; feature=1); operators, variable_names)
x2 = ComposableExpression(Node{Float64}(; feature=2); operators, variable_names)
x3 = ComposableExpression(Node{Float64}(; feature=3); operators, variable_names)
Note that using x1 here refers to the relative argument to the expression. So the node with feature equal to 1 will reference the first argument, regardless of what it is.
st_expr = TemplateExpression(
(; f=x1 - x2 * x2, g=1.5 * x1);
structure,
operators,
variable_names
) # Prints as: f = #1 - (#2 * #2); g = 1.5 * #1
# Evaluation combines evaluation of `f` and `g`, and combines them
# with the structure function:
st_expr([0.0; 1.0; 2.0;;])
This also work with hierarchical expressions! For example,
structure = TemplateStructure{(:f, :g)}(
((; f, g), (x1, x2, x3)) -> f(x1, g(x2), x3^2) - g(x3)
)
this is a valid structure!
We can also use this TemplateExpression in SymbolicRegression.jl searches!
For example, say that we want to fit *vector expressions*:
using SymbolicRegression
using MLJBase: machine, fit!, report
We first define our structure. This also has our variable mapping, which says we are fitting f(x1, x2), g1(x3), and g2(x3):
function my_structure((; f, g1, g2), (x1, x2, x3))
_f = f(x1, x2)
_g1 = g1(x3)
_g2 = g2(x3)
# We use `.x` to get the underlying vector
out = map((fi, g1i, g2i) -> (fi + g1i, fi + g2i), _f.x, _g1.x, _g2.x)
# And `.valid` to see whether the evaluations
return ValidVector(out, _f.valid && _g1.valid && _g2.valid)
end
structure = TemplateStructure{(:f, :g1, :g2)}(my_structure)
Now, our dataset is a regular 2D array of inputs for X. But our y is actually a vector of 2-tuples!
X = rand(100, 3) .* 10
y = [
(sin(X[i, 1]) + X[i, 3]^2, sin(X[i, 1]) + X[i, 3])
for i in eachindex(axes(X, 1))
]
Now, since this is a vector-valued expression, we need to specify a custom elementwise_loss function:
elementwise_loss = ((x1, x2), (y1, y2)) -> (y1 - x1)^2 + (y2 - x2)^2
This reduces y and the predicted value of y returned by the structure function.
Our regressor is then:
model = SRRegressor(;
binary_operators=(+, *),
unary_operators=(sin,),
maxsize=15,
elementwise_loss=elementwise_loss,
expression_type=TemplateExpression,
# Note - this is where we pass custom options to the expression type:
expression_options=(; structure),
)
mach = machine(model, X, y)
fit!(mach)
Let’s see the performance of the model:
report(mach)
We can also check the expression is split up correctly:
r = report(mach)
idx = r.best_idx
best_expr = r.equations[idx]
best_f = get_contents(best_expr).f
best_g1 = get_contents(best_expr).g1
best_g2 = get_contents(best_expr).g2
3. Created “Parametric Expressions”, for custom functional forms with per-class parameters: (ParametricExpression <: AbstractExpression)
Parametric Expressions are another example of an AbstractExpression with additional features than a normal Expression. This type allows SymbolicRegression.jl to fit a parametric functional form, rather than an expression with fixed constants. This is particularly useful when modeling multiple systems or categories where each may have unique parameters but share a common functional form and certain constants.
A parametric expression is constructed with a tree represented as a ParametricNode <: AbstractExpressionNode – this is an alternative type to the usual Node type which includes extra fields: is_parameter::Bool, and parameter::UInt16. A ParametricExpression wraps this type and stores the actual parameter matrix (under .metadata.parameters) as well as the parameter names (under .metadata.parameter_names).
Various internal functions have been overloaded for custom behavior when fitting parametric expressions. For example, mutate_constant will perturb a row of the parameter matrix, rather than a single parameter.
When fitting a ParametricExpression, the expression_options parameter in Options/SRRegressor should include a max_parameters keyword, which specifies the maximum number of separate parameters in the functional form.
Let's see an example of fitting a parametric expression:
using SymbolicRegression
using Random: MersenneTwister
using Zygote
using MLJBase: machine, fit!, predict, report
Let’s generate two classes of model and try to find it:
rng = MersenneTwister(0)
X = NamedTuple{(:x1, :x2, :x3, :x4, :x5)}(ntuple(_ -> randn(rng, Float32, 30), Val(5)))
X = (; X..., classes=rand(rng, 1:2, 30)) # Add class labels (1 or 2)
# Define per-class parameters
p1 = [0.0f0, 3.2f0]
p2 = [1.5f0, 0.5f0]
# Generate target variable y with class-specific parameters
y = [
2 * cos(X.x4[i] + p1[X.classes[i]]) + X.x1[i]^2 - p2[X.classes[i]]
for i in eachindex(X.classes)
]
When fitting a ParametricExpression, it tends to be more important to set up an autodiff_backend such as :Zygote or :Enzyme, as the default backend (finite differences) can be too slow for the high-dimensional parameter spaces.
model = SRRegressor(
niterations=100,
binary_operators=[+, *, /, -],
unary_operators=[cos, exp],
populations=30,
expression_type=ParametricExpression,
expression_options=(; max_parameters=2), # Allow up to 2 parameters
autodiff_backend=:Zygote, # Use Zygote for automatic differentiation
parallelism=:multithreading,
)
mach = machine(model, X, y)
fit!(mach)
The expressions are returned with the parameters:
r = report(mach);
best_expr = r.equations[r.best_idx]
@show best_expr
@show get_metadata(best_expr).parameters
4. Introduced a variety of new abstractions for user extensibility
v1 introduces several new abstract types to improve extensibility. These allow you to define custom behaviors by leveraging Julia’s multiple dispatch system.
Expression types: AbstractExpression: As explained above, SymbolicRegression now works on Expression rather than Node, by default. Actually, most internal functions in SymbolicRegression.jl are now defined on AbstractExpression, which allows overloading behavior. The expression type used can be modified with the expression_type and node_type options in Options.
expression_type: By default, this isExpression, which simply stores a binary tree (Node) as well as thevariable_names::Vector{String}andoperators::DynamicExpressions.OperatorEnum. See the implementation ofTemplateExpressionandParametricExpressionfor examples of what needs to be overloaded.node_type: By default, this will beDynamicExpressions.default_node_type(expression_type), which allowsParametricExpressionto default toParametricNodeas the underlying node type.
Mutation types: mutate!(tree::N, member::P, ::Val{S}, mutation_weights::AbstractMutationWeights, options::AbstractOptions; kws...) where {N<:AbstractExpression,P<:PopMember,S}, where S is a symbol representing the type of mutation to perform (where the symbols are taken from the mutation_weights fields). This allows you to define new mutation types by subtyping AbstractMutationWeights and creating new mutate! methods (simply pass the mutation_weights instance to Options or SRRegressor).
Search states: AbstractSearchState: this is the abstract type for SearchState which stores the search process’s state (such as the populations and halls of fame). For advanced users, you may wish to overload this to store additional state details. (For example, LaSR stores some history of the search process to feed the language model.)
Global options and full customization: AbstractOptions and AbstractRuntimeOptions: Many functions throughout SymbolicRegression.jl take AbstractOptions as an input. The default assumed implementation is Options. However, you can subtype AbstractOptions to overload certain behavior.
For example, if we have new options that we want to add to Options:
Base.@kwdef struct MyNewOptions
a::Float64 = 1.0
b::Int = 3
end
we can create a combined options type that forwards properties to each corresponding type:
struct MyOptions{O<:SymbolicRegression.Options} <: SymbolicRegression.AbstractOptions
new_options::MyNewOptions
sr_options::O
end
const NEW_OPTIONS_KEYS = fieldnames(MyNewOptions)
# Constructor with both sets of parameters:
function MyOptions(; kws...)
new_options_keys = filter(k -> k in NEW_OPTIONS_KEYS, keys(kws))
new_options = MyNewOptions(; NamedTuple(new_options_keys .=> Tuple(kws[k] for k in new_options_keys))...)
sr_options_keys = filter(k -> !(k in NEW_OPTIONS_KEYS), keys(kws))
sr_options = SymbolicRegression.Options(; NamedTuple(sr_options_keys .=> Tuple(kws[k] for k in sr_options_keys))...)
return MyOptions(new_options, sr_options)
end
# Make all `Options` available while also making `new_options` accessible
function Base.getproperty(options::MyOptions, k::Symbol)
if k in NEW_OPTIONS_KEYS
return getproperty(getfield(options, :new_options), k)
else
return getproperty(getfield(options, :sr_options), k)
end
end
Base.propertynames(options::MyOptions) = (NEW_OPTIONS_KEYS..., fieldnames(SymbolicRegression.Options)...)
These new abstractions provide users with greater flexibility in defining the structure and behavior of expressions, nodes, and the search process itself. These are also of course used as the basis for alternate behavior such as ParametricExpression and TemplateExpression.
5. Added TensorBoardLogger.jl and other logging integrations via SRLogger
You can now track the progress of symbolic regression searches using TensorBoardLogger.jl, Wandb.jl, or other logging backends.
This is done by wrapping any AbstractLogger with the new SRLogger type, and passing it to the logger option in SRRegressor or equation_search:
using SymbolicRegression
using TensorBoardLogger
logger = SRLogger(
TBLogger("logs/run"),
log_interval=2, # Log every 2 steps
)
model = SRRegressor(;
binary_operators=[+, -, *],
logger=logger,
)
The logger will track:
- Loss curves over time at each complexity level
- Population statistics (distribution of complexities)
- Pareto frontier volume (can be used as an overall metric of search performance)
- Full equations at each complexity level
This works with any logger that implements the Julia logging interface.
6. Support for Zygote.jl and Enzyme.jl within the constant optimizer, specified using the autodiff_backend option
Historically, SymbolicRegression has mostly relied on finite differences to estimate derivatives – which actually works well for small numbers of parameters. This is what Optim.jl selects unless you can provide it with gradients.
However, with the introduction of ParametricExpressions, full support for autodiff-within-Optim.jl was needed. v1 includes support for some parts of DifferentiationInterface.jl, allowing you to actually turn on various automatic differentiation backends when optimizing constants. For example, you can use
Options(
autodiff_backend=:Zygote,
)
to use Zygote.jl for autodiff during BFGS optimization, or even
Options(
autodiff_backend=:Enzyme,
)
for Enzyme.jl (though Enzyme support is highly experimental).
7. Other Changes
Changed output file handling
Instead of writing to a single file like hall_of_fame_<timestamp>.csv, outputs are now organized in a directory structure. Each run gets a unique ID (containing a timestamp and random string, e.g., 20240315_120000_x7k92p), and outputs are saved to outputs/<run_id>/. Currently, only saves hall_of_fame.csv (and hall_of_fame.csv.bak), with plans to add more logs and diagnostics in this folder in future releases.
The output directory can be customized via the output_directory option (defaults to ./outputs). A custom run ID can be specified via the new run_id parameter passed to equation_search (or SRRegressor).
Other Small Features in v1.0.0
- Support for per-variable complexity, via the
complexity_of_variablesoption. - Option to force dimensionless constants when fitting with dimensional constraints, via the
dimensionless_constants_onlyoption. - Default
maxsizeincreased from 20 to 30. - Default
niterationsincreased from 10 to 50, as many users seem to be unaware that this is small (and meant for testing), even in publications. I think this 50 is still low, but it should be a more accurate default for those who don’t tune. MLJ.fit!(mach)now records the number of iterations used, and, shouldmach.model.niterationsbe changed after the fit, the number of iterations passed toequation_searchwill be reduced accordingly.- Fundamental improvements to the underlying evolutionary algorithm
- New mutation operators introduced,
swap_operandsandrotate_tree– both of which seem to help kick the evolution out of local optima. - New hyperparameter defaults created, based on a Pareto front volume calculation, rather than simply accuracy of the best expression.
- New mutation operators introduced,
- Major refactoring of the codebase to improve readability and modularity
- Identified and fixed a major internal bug involving unexpected aliasing produced by the crossover operator
- Segmentation faults caused by this are a likely culprit for some crashes reported during multi-day multi-node searches.
- Introduced a new test for aliasing throughout the entire search state to prevent this from happening again.
- Improved progress bar and StyledStrings integration.
- Julia 1.10 is now the minimum supported Julia version.
- Also see the “Update Guide” below for more details on upgrading.
Update Guide
Note that most code should work without changes! Only if you are interacting with the return types of equation_search or report(mach), or if you have modified any internals, should you need to make some changes.
Also note that the “hall of fame” CSV file is now stored in a directory structure, of the form outputs/<run_id>/hall_of_fame.csv. This is to accommodate additional log files without polluting the current working directory. Multi-output runs are now stored in the format .../hall_of_fame_output1.csv, rather than the old format hall_of_fame_{timestamp}.csv.out1.
So, the key changes are, as discussed above, the change from Node to Expression as the default type for representing expressions. This includes the hall of fame object returned by equation_search, as well as the vector of expressions stored in report(mach).equations for the MLJ interface. If you need to interact with the internal tree structure, you can use get_contents(expression) (which returns the tree of an Expression, or the named tuple of a ParametricExpression - use get_tree to map it to a single tree format).
To access other info stored in expressions, such as the operators or variable names, use get_metadata(expression).
This also means that expressions are now basically self-contained. Functions like eval_tree_array no longer require options as arguments (though you can pass it to override the expression’s stored options). This means you can also simply call the expression directly with input data (in [n_features, n_rows] format).
Before this change, you might have written something like this:
using SymbolicRegression
x1 = Node{Float64}(; feature=1)
options = Options(; binary_operators=(+, *))
tree = x1 * x1
This had worked, but only because of some spooky action at a distance behavior involving a global store of last-used operators! (Noting that Node simply stores an index to the operator to be lightweight.)
After this change, things are much cleaner:
options = Options(; binary_operators=(+, *))
operators = options.operators
variable_names = ["x1"]
x1 = Expression(Node{Float64}(; feature=1); operators, variable_names)
tree = x1 * x1
This is now a safe and explicit construction, since * can lookup what operators each expression uses, and infer the right indices! This operators::OperatorEnum is a tuple of functions, so does not incur dispatch costs at runtime. (The variable_names is optional, and gets stripped during the evolution process, but is embedded when returned to the user.)
We can now use this directly:
println(tree) # Uses the `variable_names`, if stored
tree(randn(1, 50)) # Evaluates the expression using the stored operators
Also note that the minimum supported version of Julia is now 1.10. This is because Julia 1.9 and earlier have now reached end-of-life status, and 1.10 is the new LTS release.
Additional Notes
- Custom Loss Functions: Continue to define these on
AbstractExpressionNode. - General Usage: Most existing code should work with minimal changes.
- CI Updates: Tests are now split into parts for faster runs, and use TestItems.jl for better scoping of test variables.
Be sure to check out the documentation here: Home · SymbolicRegression.jl
You can also see the full CHANGELOG.md for a history of the library.
In general I’m very excited to see what people can do with these new features. They make SymbolicRegression.jl a lot more modular and extensible, so heavy customizations are now possible.