GP_ELITE is a symbolic regression engine in pure Python: given (X, y) data, it
searches for a readable mathematical formula linking them, instead of a
black-box model.
To show what that means concretely: I gave it nothing but the 8 planets'
distance from the Sun and orbital period — 8 data points — and asked for a
formula. It returned:
T = a · sqrt(a) (i.e. a^1.5), R² = 1.000000
That's Kepler's Third Law (T² ∝ a³), which took Kepler ~10 years to find in
1618. GP_ELITE found it in ~3 seconds. Reproducible: examples/kepler_demo.py.
v0.2.0 (this week) added the parts that make it reliable: Levenberg-Marquardt
constant fitting (constants come back at machine precision — Coulomb's
q1·q2/(4πεr²) is recovered exactly), multi-restart with a merged candidate
archive, a Pareto front output (the full complexity ↔ accuracy staircase, not
just one champion), and a guarded forecasting mode for extrapolating trends
beyond your data without the usual GP blow-ups.
Pure Python/NumPy — pip install gp-elite, no compiler, no Julia.
Target Audience
Anyone with small experimental datasets (≤10 variables, 100–5000 points) who
wants to understand a relationship, not just predict it: lab engineers,
scientists, students. One concrete use case that drove development: battery
degradation (SOH) forecasting — the guarded mode gives you an honest bracket
of scenarios (a Pareto front from a conservative straight line to richer
bounded laws) instead of one overconfident curve. Production-usable for that
niche (built-in hold-out validation, regression-tested); not aimed at
large-scale ML.
Comparison
vs gplearn (the established pure-Python option): I ran both on the same
frozen benchmark — 15 Feynman physics equations, identical data and splits,
generous budget for gplearn. Exact symbolic recovery (machine precision):
GP_ELITE 10/15 (67%) vs gplearn 6/15 (40%). gplearn recovers the
constant-free formulas and stalls as soon as a ½ or a 4π appears (no real
constant optimization); LM fitting is what closes that gap. Every number is
reproducible: PYTHONHASHSEED=0 python benchmarks/feynman_bench.py 0 15 and
benchmarks/duel.py in the repo.
vs PySR / Operon (the state of the art): they are stronger on speed and
scale, and I'm not claiming otherwise — but they require a Julia or C++
toolchain. GP_ELITE's whole point is zero barrier: pip install and go.
vs neural nets / gradient boosting: those win on raw accuracy for large
data, but give you a black box — GP_ELITE gives you the actual equation.
Honest limits: weak on chaotic targets (tested on Collatz), degrades past ~6
variables with decoy features, and pure Python costs wall-time on big data.
This is not surprising at all and depends on the inductive bias hardcoded in the search.
There are infinite number of curves that agree on those 8 points and deviate from Kepler 's law everywhere else. On such 'trajectories' this algorithm would have performed badly.
Good questions — I ran all three:
Semi-minor axis: I added b as a second feature (b = a·sqrt(1−e²), so corr(a,b) = 1.00000 to 5 decimals; Mercury is the only planet where they differ by more than 2%). It still picked a and ignored b completely: T = 164.78·a·sqrt(a), R² = 0.99999998. What saves it on such a nasty collinear decoy is the constant fitting: the law is exact in a and only almost-exact in b, so Levenberg-Marquardt makes that 2% Mercury error decisive.
Kepler's 6 planets: works fine, and it's actually nicer — it returned pow(a, 1.500812) explicitly. Six clean points on a power law is plenty.
Noise: this is where I have to be honest. With 1% gaussian noise on T the fit is still R² = 0.9996 and a·sqrt(a) is still in there, but wrapped in junk (exp(tanh(...))). At 5% the clean form is gone — it returns a bounded exp(−a²) family that fits well but isn't the law. So fit quality degrades gracefully, symbolic recovery doesn't. Making that part noise-robust is pretty much the open frontier of the whole field, not just of my tool.
> b = a·sqrt(1−e²), so corr(a,b) = 1.00000 to 5 decimals
Isn't e different for each planet?
> Mercury is the only planet where they differ by more than 2%
I remember something about Mars been the planet with the most eccentric elipse
> *So fit quality degrades gracefully, symbolic recovery doesn't. Making that part noise-robust is pretty much the open frontier of the whole field, not just of my tool.
Nice. It's a hard problem. Which heuristic are you using to pick the "best" formula?
Fair point, and nyrikki has it right: the engine never sees units. Everything is normalized to dimensionless magnitudes before the search, so that formula is an empirical fit on pure numbers, not a physical law. I should say that more clearly in the README.
Enforcing dimensional consistency during the search is a known extension (AI Feynman does a version of it) and honestly it's a good roadmap item, a dimensionally sound form would be more trustworthy.
I've seen crazy stuff for heat exchangers and other stuff used in factories.
When the temperature difference is small, everithing is linear and you get nice formulas.
When the difference of temperature is big and you have liquids with convection and turbulence, you only get empirical formulas with weird exponents. The correct method would be to give the constants with the correct units, but it's usual to specify unit to measure the data instead, and just enter the numbers in the formula.
And I wonder if somebody has tried with the available galactic data and see if the genetic programming can come up with a better formula than MOND or Einstein's general relativity.
For simple problems as Kepler's law, a quick detour on Desmos will show a perfect fit for power law instantly. In general, there are many important criteria for a better curve fitting (for ex. independent, normal distributed residuals), not just R, so I hope the author has/will incorporate them into the search to create a more robust result.
Funny you ask, I actually spent a weekend on exactly that. Took the SPARC rotation curve data (~2700 points, 150+ galaxies), built the radial acceleration relation from the raw files and let the engine search. Honest result: it recovered a0 around 1.2e-10 m/s2 like the published fits, and found a different functional form (a log-parabola) that fits exactly as well as the MOND interpolating function, difference of 0.0004 dex. But "exactly as well" is the problem: at the current scatter the data can't tell the forms apart. Which I guess is why the interpolating function debate never dies. No new physics from me, but it was a good way to hit the real wall of the field.
And agreed on residuals, R2 alone is weak. Hold-out validation is already in there, proper residual diagnostics are a fair ask for the roadmap.
13 comments
[ 1.2 ms ] story [ 47.8 ms ] threadGP_ELITE is a symbolic regression engine in pure Python: given (X, y) data, it searches for a readable mathematical formula linking them, instead of a black-box model.
To show what that means concretely: I gave it nothing but the 8 planets' distance from the Sun and orbital period — 8 data points — and asked for a formula. It returned:
T = a · sqrt(a) (i.e. a^1.5), R² = 1.000000
That's Kepler's Third Law (T² ∝ a³), which took Kepler ~10 years to find in 1618. GP_ELITE found it in ~3 seconds. Reproducible: examples/kepler_demo.py.
v0.2.0 (this week) added the parts that make it reliable: Levenberg-Marquardt constant fitting (constants come back at machine precision — Coulomb's q1·q2/(4πεr²) is recovered exactly), multi-restart with a merged candidate archive, a Pareto front output (the full complexity ↔ accuracy staircase, not just one champion), and a guarded forecasting mode for extrapolating trends beyond your data without the usual GP blow-ups.
Pure Python/NumPy — pip install gp-elite, no compiler, no Julia.
Target Audience
Anyone with small experimental datasets (≤10 variables, 100–5000 points) who wants to understand a relationship, not just predict it: lab engineers, scientists, students. One concrete use case that drove development: battery degradation (SOH) forecasting — the guarded mode gives you an honest bracket of scenarios (a Pareto front from a conservative straight line to richer bounded laws) instead of one overconfident curve. Production-usable for that niche (built-in hold-out validation, regression-tested); not aimed at large-scale ML.
Comparison
vs gplearn (the established pure-Python option): I ran both on the same frozen benchmark — 15 Feynman physics equations, identical data and splits, generous budget for gplearn. Exact symbolic recovery (machine precision): GP_ELITE 10/15 (67%) vs gplearn 6/15 (40%). gplearn recovers the constant-free formulas and stalls as soon as a ½ or a 4π appears (no real constant optimization); LM fitting is what closes that gap. Every number is reproducible: PYTHONHASHSEED=0 python benchmarks/feynman_bench.py 0 15 and benchmarks/duel.py in the repo.
vs PySR / Operon (the state of the art): they are stronger on speed and scale, and I'm not claiming otherwise — but they require a Julia or C++ toolchain. GP_ELITE's whole point is zero barrier: pip install and go.
vs neural nets / gradient boosting: those win on raw accuracy for large data, but give you a black box — GP_ELITE gives you the actual equation.
Honest limits: weak on chaotic targets (tested on Collatz), degrades past ~6 variables with decoy features, and pure Python costs wall-time on big data.
Code (MIT): https://github.com/ariel95500-create/gp-elite
There are infinite number of curves that agree on those 8 points and deviate from Kepler 's law everywhere else. On such 'trajectories' this algorithm would have performed badly.
What happens if you give the system not only the semi-mayor axis but also the semi-minor axis?
Have you tried with only the 6 planets Kepler know? (I don't expect this to change the result too much.)
Have you tired with noisy data?
to get a new paragraph here.
> b = a·sqrt(1−e²), so corr(a,b) = 1.00000 to 5 decimals
Isn't e different for each planet?
> Mercury is the only planet where they differ by more than 2%
I remember something about Mars been the planet with the most eccentric elipse
> *So fit quality degrades gracefully, symbolic recovery doesn't. Making that part noise-robust is pretty much the open frontier of the whole field, not just of my tool.
Nice. It's a hard problem. Which heuristic are you using to pick the "best" formula?
capacity_SOH ≈ 0.913 − 0.352 · tanh( cycle^((temperature/cycle)^0.485) )
I understand this fits the data, but exponents should be dimensionless, what is temperature/cycle?
Enforcing dimensional consistency during the search is a known extension (AI Feynman does a version of it) and honestly it's a good roadmap item, a dimensionally sound form would be more trustworthy.
When the temperature difference is small, everithing is linear and you get nice formulas.
When the difference of temperature is big and you have liquids with convection and turbulence, you only get empirical formulas with weird exponents. The correct method would be to give the constants with the correct units, but it's usual to specify unit to measure the data instead, and just enter the numbers in the formula.
After a short search in Google, I got this example https://www.researchgate.net/figure/Equations-used-to-obtain...
For simple problems as Kepler's law, a quick detour on Desmos will show a perfect fit for power law instantly. In general, there are many important criteria for a better curve fitting (for ex. independent, normal distributed residuals), not just R, so I hope the author has/will incorporate them into the search to create a more robust result.
And agreed on residuals, R2 alone is weak. Hold-out validation is already in there, proper residual diagnostics are a fair ask for the roadmap.