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

25 pointsposted 6 days ago
by sade_95
(github.com)

11 Comments

sinuhe69

13 minutes ago

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.

ziofill

an hour ago

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?

nyrikki

11 minutes ago

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...

gus_massa

4 days ago

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?

sade_95

2 days ago

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.

gus_massa

2 days ago

Remember to use two enters

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?

sade_95

2 days ago

Thanks for the formatting tip, noted.

Isn't e different for each planet?

Yes, each planet got its own e (Mercury 0.206, Venus 0.007, Earth 0.017...). The correlation still comes out at 1.00000 because all eccentricities are small while a spans 0.39 to 30 AU — a ≤2% per-planet wobble is invisible to Pearson across two decades of range. That's exactly what makes it a nasty decoy: almost collinear, but not quite.

Mars been the planet with the most eccentric ellipse

Close — Mercury is actually the most eccentric (0.206), Mars is second (0.093). Funny enough, Mars is the famous one precisely because Kepler derived his laws fighting with Tycho's Mars data: its ellipse was just eccentric enough to kill every circular fit. If Tycho had handed him Venus data instead, we might have waited a while longer.

Which heuristic are you using to pick the "best" formula?

Hold-out validation during evolution, then the old CART-style one-standard-error rule: pick the smallest formula whose validation error is within ~1 SE of the best one. Constants get refit with Levenberg-Marquardt before that comparison, so small forms compete at their best. It also returns the full accuracy-vs-size Pareto front so you can override the choice. And to connect it to your noise question: under noise that 1-SE band is exactly where wrong-but-simpler formulas sneak in and tie the true one — the two problems are really the same problem.

sade_95

6 days ago

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

srean

2 hours ago

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.

sade_95

2 hours ago

You're right, and I'd go further: the inductive bias isn't incidental, it's the whole product. Short trees over {+, *, sqrt, exp, ...} plus a parsimony penalty is basically Occam's razor made executable. A bias-free learner can't generalize at all (no free lunch), so the honest question is whether this particular bias matches the domain. For physics it has an unreasonably good track record though that's the mystery of physics, not of my library. Two nuances though. The evidence here isn't "a curve fits 8 points" infinitely many do, as you say. It's that a 3-node formula fits them to machine precision (1−R² ≈ 1e-15). Under an MDL view that's not nothing: the probability that such a short description nails 8 independent points exactly, if the truth were some unrelated wiggly curve, is astronomically small. The shortness is the evidence. Second: your adversarial curve would fool Kepler too, and any finite-data method ever. The practical mitigation is the boring one held-out validation, and in the planetary case, extrapolation: the law found on 8 planets keeps working on moons, asteroids and exoplanets. On the tool side I try to keep the failure mode visible rather than hidden: it returns the full accuracy-vs-complexity Pareto front, and the docs say plainly that noise breaks symbolic recovery long before it breaks fit quality. So yes: it finds simple laws when simple laws exist. When they don't, it fails ideally loudly.