comparion between python and julia version of symbolic regression

[wenpw]

Dear Miles,

Thank you very much for the powerful code.
Please help check the following problems.

Recently, I have compared the Pysr and the Julia version to fit the same formula:
y = 4log(1+exp(2pi*(x-100)/4)/2/x.

The codes are in the attached zip file.

Two serious problems were found:

In Julia, it takes only few minutes to obtain the correct formula using only two cpus.
In Python, I have tried for more than one day using even twenty cpu, however, I still could not get the right result.
I have also tried different parameters and could still not find the correct answer.

May I ask why this happens, and how to improve the python script to get the correct results?

In julia, hundreds of G RAM memory are needed in my servers for the running for example using numprocs=20.
In python, no large RAM memroy is needed in my test.
Is there any way to release the memory during the running of Julia code?

Thank you very much in advance.

log1p-v2.zip

[MilesCranmer]

Dear @wenpw,
Thanks for posting. It looks like your search settings are different in Julia vs Python. Could you please try with the same settings?

Also can you paste the code into your answer directly rather than through a zip file?
Thanks,
Miles

[wenpw]

Dear @MilesCranmer:

Thank you for the quick reply.
In fact, I have tried many different settings for these two version.
And the problems are smilar.

Here are the codes:

1. The python version:

import numpy as np
import sympy

def wong(expe,vb,rb,curv):
    fwong =  curv  * np.log(1.0 +np.exp(2.*np.pi*(expe-vb)/curv)) /2./expe
    return fwong

x0 = 85+np.random.rand(300,1)*150.
x1 = 100.
x2 = 10.
x3 = 4.

curv = x3*np.ones((300,1))
rb = x2*np.ones((300,1))
vb = x1*np.ones((300,1))
expe = x0
y = wong(expe,vb,rb,curv)

X = expe

from pysr import PySRRegressor

model = PySRRegressor(
    procs=400,
    populations=200,
    # ^ 2 populations per core, so one is always running.
    population_size=500,
    # ncyclesperiteration=550, 
    # ^ Generations between migrations.
    niterations=10000000,  # Run forever
    early_stop_condition=(
        "stop_if(loss, complexity) = loss < 1e-6 && complexity < 10"
        # Stop early if we find a good and simple equation
    ),
    timeout_in_seconds=60 * 60 * 24 * 20,
    maxsize=20,
    binary_operators=["*", "+", "-", "/"],
    unary_operators=[
        "exp",
        "log1p",
    ],
    loss="loss(x, y) = abs(x - y)^2 + 30000*(log(abs(x)) - log(abs(y)))^2",
    # loss="loss(x, y) =  (log(abs(x)+1e-48) - log(abs(y)+1e-60))^2",
    # ^ Custom loss function (julia syntax)
    progress=True,
    warm_start=True,
    turbo=False,
)

model.fit(X, y)

print(model)

print(model.predict(X))

2. the Julia version:

using SymbolicRegression, SymbolicUtils

expe = 85.0f0 .+ rand(Float64, 1, 300) .*150#.0f0
curv = expe[1,:] .- expe[1,:] .+ 4.0f0
rb = expe[1,:] .- expe[1,:] .+ 10.0f0
vb = expe[1,:] .- expe[1,:] .+ 100.0f0


X=expe
y =  curv[1,:] .*  log.(1.0f0 .+ exp.(2.0f0 * pi .*(expe[1,:]  .-vb[1,:] )./ curv[1,:])  ) ./2.0f0 ./expe[1,:]


early_stop(loss, c) = ((loss <= 1e-10) && (c <= 10))
my_loss(x, y) = abs(x - y)^2 + 30000*(log(abs(x)) - log(abs(y)))^2

options = SymbolicRegression.Options(;
    binary_operators=[+, *, /, -], 
    unary_operators=[exp, log1p], 
    npopulations=200,
    early_stop_condition=early_stop,
    enable_autodiff=true,
    elementwise_loss=my_loss,
)

hall_of_fame = EquationSearch(
    X, y;
    niterations=400000000, 
    options=options, 
    numprocs=2,
    parallelism=:multiprocessing,
)

dominating = calculate_pareto_frontier(hall_of_fame)

trees = [member.tree for member in dominating]

tree = trees[end]
output, did_succeed = eval_tree_array(tree, X, options)

eqn = node_to_symbolic(dominating[end].tree, options)
println("Complexity\tMSE\tEquation")

for member in dominating
    complexity = compute_complexity(member.tree, options)
    loss = member.loss
    string = string_tree(member.tree, options)
    println("$(complexity)\t$(loss)\t$(string)")
end

[MilesCranmer]

(e.g., you requested 400 cores in PySR, but only 2 cores in SymbolicRegression.jl)

[wenpw]

Year, the first serious problem is:

Even when I use 400 cores in Pysr, I could not get the correct answer for more than 1 day's running.
Using 2 cores, the answer could still not be obtained.

In julia, 2 cores are encough within few minuits but with huge RAM memories.

[MilesCranmer]

Even when I use 400 cores in Pysr, I could not get the correct answer for more than 1 day's running.

Just to be clear – if you only have 20 cores, then you should not set a higher number than 20 for procs. Even if it is working, that would slow things down quite a bit as it overload the CPU scheduler.

[MilesCranmer]

I noticed you have parallelism=:multiprocessing in the Julia version. Could you try using that for Python as well? In the Python version, that would be multithreading=False.

[MilesCranmer]

Also note that the way you use 64-bit precision in PySR is by passing precision=64, not by passing data of 64-bit type. This is one difference between PySR and SymbolicRegression.jl's frontends.

[wenpw]

Thank you for your patience.

Can you share details on your measurements of execution speed? For example, a screenshot of htop? The Ubuntu data will be more useful for me.

Here is the htop data in my Ubuntu server:

For the record I'm not able to reproduce your measurements; your PySR example runs fast for me on both macOS and linux:

Do you mean that you could obtain the corrrect answer by Pysr on macOS and Linux very quickly?

[wenpw]

Also note that the way you use 64-bit precision in PySR is by passing precision=64, not by passing data of 64-bit type. This is one difference between PySR and SymbolicRegression.jl's frontends.

Thank you for this reminding.

[MilesCranmer]

Wait I am confused. I thought you mentioned you were not able to get multithreading working? It looks like it is using multiple cores there, no?

[wenpw]

Yeah, there is no problems using multiple cores. The Problem is using Pysr, no matter how I adjust the parameters, I could not get the correct answer. Using julia, correct answer could be obtained very soon even with only 1 cpu.