Pure-Python symbolic regression that rediscovered Kepler's law from 8 data point

Join the discussion

• Hacker News
• What My Project Does

  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.

  Code (MIT): https://github.com/ariel95500-create/gp-elite

  by sade_95
• 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.

  by srean
• Total slop from the first word
  by etdznots
• 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.

  by sinuhe69
• Some minor comments:

  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?

  by gus_massa
• 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.
  by sade_95
• The battery example makes no sense:

  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?

  by ziofill
• I am not saying this result is correct, but you don't have unit safety in python at all

  In this example temperature would be a magnitude, not a unitful value.

  At least in the ISO 8000whatever convention where a value is the product of a unit and a magnitude like most people are use to.

  Here is a Terry Tao post with more information[0] on why the convention is there, but as he mentions, in differential geometry and Clifford/Geometric Algebra you do things like add vectors to scalers all the time.

  [0] https://terrytao.wordpress.com/2012/12/29/a-mathematical-for...

  by nyrikki