A while ago I wrote The Hidden Geometry of Numbers: every whole number is a point in a lattice with one axis per prime, multiplying is adding arrows, and a sheet of glass sliding through that lattice finds the primes without dividing anything.

I always thought that it would be such a nice experience to watch it in a video, but making one takes a lot of effort. In this article I tell you what tools I used and how I solved the main problems.

What’s in the video

The video opens with a story: around 500 BC, in what is today Crotone in southern Italy, Pythagoras gathered a brotherhood around one belief: everything is number, and every length is a ratio of whole numbers. One of them, Hippasus, looked at the diagonal of a square with sides of length one is not a ratio of whole numbers. To a brotherhood whose motto was all is number, that was heresy. Legend says he went out to sea and never came back.

From there the video follows the article: numbers as places on a grid, multiplication as adding arrows, divisors, gcd and lcm, a new dimension for every prime, and fractions filling the other quadrants. Halfway through it returns to Hippasus, because by then his proof can be seen at a glance. The square root of 2 sits half a step along the 2-axis, between the grid points, while every fraction sits on a grid point. The last part is the ending of the article: the sheet of glass that sweeps through the lattice, catches the primes, and touches exactly one point for every number.

The hardest part of an explainer video is timing

When the narrator says “lay six’s arrow on its tip”, the arrow has to move at that moment, not half a second later. Lining this up by hand for fourteen minutes would be a lot of work, and every change to a sentence would shift everything after it.

To prevent myself from the timing hell above, I used the following approach: the voice was made first, the timing was taken from the voice, and the animation was written last, to follow that timing. This way I could avoid figuring out timings of an animation and just time it to the audio generated, even when I made changes in the narration. And I solved it in Python.

The script

The script is a Python file rather than a document: 14 scenes, 102 lines and 1,979 spoken words. Wherever something should happen on the screen, the line contains a mark in curly braces. Shortened, one of them looks like this:

f"Take seventy-two's arrow, {g6}and lay six's arrow on its tip."

The mark {g6} stands for the moment the narrator starts the next word, in this case "and". The whole script has 173 marks of this kind; a few of them only ask for a pause.

The voice

Each line was sent separately to a voice model, called Fish Audio Drama 3, running on my computer. A speech recognizer, Whisper, then listened to every recording and compared it word by word with the script. When a word was swallowed or missing, the line was recorded again. The final narration uses 148 recordings for the 102 lines, because 28 lines needed more than one attempt. Generating and checking the recordings that made it into the video took about 70 minutes.

The timing

Whisper gives more than the words it heard: it also says when each word starts and ends. So once the voice existed, I knew to a hundredth of a second when the narrator says the “and” in the sentence above: 62.13 seconds into its scene. The recordings were then joined into one narration track, with the pauses the script asks for, and a timing file was written next to it. That file records when each scene starts, when each line starts, and when each mark happens.

The animation

Only then came the picture. The animations are made with Manim, the library 3Blue1Brown built for his videos, in its Community edition. The scene code never refers to a time in seconds. It waits for a mark and then moves. Simplified, the moment above reads:

self.at("g6")    # wait until the narrator says "and"
self.play(...)   # then lay six's arrow on seventy-two's tip

During rendering, the scene looks up g6 in the timing file and waits exactly that long. This is what keeps the code independent of any particular recording. When I rewrote a sentence, only that line had to be recorded again; the timing file was rebuilt, and every animation after it moved to its new place without any manual adjustment.

One technical detail: Manim has a proper 3-D mode, but it was far too slow for clouds of glowing points. In a test it took about five minutes to render two seconds of film. The 3-D scenes are therefore ordinary 2-D scenes with a small perspective camera written by hand.

Sound and the final cut

Finally, the narration track was mixed with a quiet music bed that steps back whenever the narrator speaks, and with a few sound effects generated in code. The plucked string at the beginning, for example, is a few lines of Karplus–Strong synthesis. The loudness was normalized to −14 LUFS, and the rendered scenes were joined with short crossfades.

Reviewing the cuts

While the video was being made, I watched it as a series of rough cuts: low-resolution previews at 480p, with sound, which render quickly. After each one I wrote notes, a code change3, and a new cut came back. Only when no notes were left would the video be rendered once at full quality, 1080p at 60 frames per second.

Some of the most important notes were about tone rather than mathematics. The voice model accepts a short instruction for every line, describing how to read it. For the dramatic moments of the Hippasus story the first instructions asked for a hushed delivery, and the model took that literally and whispered. It sounded strange rather than dramatic, so every line went back to a plain, calm narration. The opposite problem appeared in a few sentences written with too much excitement, such as “you just multiplied two numbers without multiplying anything!”, which sounded like a children’s programme. Those sentences were rewritten in a more even tone.

The most useful note was about a picture. The first 3-D scenes looked impressive in motion, but they were hard to read. When I put a frame from the video next to the plain matplotlib figure from my own article, the article’s version was clearly easier to understand. It shows the lattice inside a room with gridded walls, tick marks and labelled points, and those give the eye something to measure depth against. The 3-D scenes of the video were rebuilt in that style: the same room, a camera with almost no perspective distortion, sharp points, and a slow rotation that helps the eye read the depth.

Problems along the way

A few problems were instructive enough to describe.

The string that stopped vibrating

The video opens with a plucked string playing an octave and a fifth, the whole-number ratios on which the Pythagoreans built their belief. The string’s movement on screen is driven by the actual loudness of the synthesized sound, so it should vibrate for exactly as long as the note can be heard. In the first cuts it didn’t: it moved briefly after each pluck and then froze while the note was still sounding.

Why did it happen? Manim tries to save work during pauses. The documentation of Scene.wait says of its frozen_frame parameter: "If None (the default), the scene tries to determine whether or not the frame is frozen on its own." It decides this by looking for updaters that depend on elapsed time. The string's movement was computed from the scene clock in a different way, so Manim concluded that nothing was moving and repeated the same frame. The fix was a single argument:

self.wait(dt, frozen_frame=False)

The “ch” that was taken for noise

Voice models sometimes leave a short breath or click after the end of a line, and at a cut this is heard as a glitch. The pipeline therefore cleaned the end of every recording: a sound that came after a silence was treated as noise and removed.

Then I noticed that the sentence “…predicting which point its sheet will touch” ended in “tou-”. The word “touch” contains a short silence of its own, the moment when the tongue closes before the “ch” is released. The cleaning step took that silence for the end of the line and removed the “ch” as noise.

The rule was then made stricter. A sound now counts as noise only if at least 0.35 seconds of silence separate it from the speech and it is shorter than 0.6 seconds. The pause inside a consonant is never that long.

A file that Windows would not play

The scenes are joined with crossfades in ffmpeg. The first finished cut would not open in the Windows Media Player, which reported “It uses unsupported encoding settings. 0x80004005.” It turned out that the crossfade filter had switched the video to the yuv444p pixel format, which many players and hardware decoders don’t support. Converting back to yuv420p at the end of the filter chain solved it. If a video made with ffmpeg plays in one player but not in another, the pixel format is a good first thing to check.

The thumbnail

Image generators are very good at atmosphere and poor at precision. A request for “a lattice with points labelled 1 to 6 and a plane passing them” produces something attractive but mathematically wrong.

Some of the thumbnails were made in two layers. ChatGPT generated only the backgrounds, from prompts that ruled out any text, kept all geometry soft and out of focus, and described which parts of the image to leave free: the top for the title, the left side for the diagram, and the right side for the red 7. Manim then drew everything that had to be correct on a transparent background: the room, the lattice, the glass sheet at the height of 7, the labels and the title. A short Python step placed the drawing over the background and added a soft glow.

Two details mattered. First, the position of the drawing was measured from the background rather than estimated: the centre of the starburst was located in the pixels, and the lattice’s origin was placed exactly on it. Second, the transparent images from Manim store their colours in premultiplied form, and the fully transparent pixels still carried the scene’s background colour. Combined in the usual way, every semi-transparent edge became too dark. Rendering the drawing over black and combining it as premultiplied colour fixed that.

Looking back

Working this way, the voice and the timing turned out to be the backbone of the whole project. Once the animation followed the narration instead of the other way round, changing the script stopped being expensive, and that made it possible to keep improving the video through many rounds of notes.

The video is on YouTube, and the written version, with all its original figures, is here.

Pick a number. Say 12. Break it into primes — 12 = 2 × 2 × 3 — and that’s the only way. This is the fundamental theorem of arithmetic: every whole number has exactly one prime recipe.

That sounds like a fact about factoring. I say it’s really about geometry — and once you see it, you can’t unsee it. Stay to the end and you’ll watch an infinite sheet slice through an infinite lattice and kiss a single point, and catch a prime in the very act of being prime:

That’s the destination. Here’s the road to it.

From a number to a coordinate

Here’s the move. Take the factorization of 12 and write it with explicit exponents, one slot per prime:

Now throw away everything except the exponents and keep them in order:

That tuple is 12. No information was lost — uniqueness (the fundamental theorem) is exactly the promise that this address belongs to 12 and to nothing else. We’ve turned a number into a coordinate: the first axis counts powers of 2, the second counts powers of 3, the third counts powers of 5, and so on, forever. There is one dimension for every prime.

So integers live in an infinite-dimensional lattice, with the primes as the axes. That sounds intimidating, but almost every number you’d want to draw only lights up a few axes. If a number’s only prime factors are 2 and 3, all its other coordinates are zero — so it lives entirely in the 2–3 plane, and we can draw it on ordinary graph paper.

Watch the small numbers take their seats. Horizontal axis = how many 2s, vertical axis = how many 3s:

Read off a few to convince yourself:

• 1 sits at the origin (0,0) — no primes at all, the empty product.
• 2 is (1,0), 4 is (2,0), 8 is (3,0): walking right doubles you each step.
• 3 is (0,1), 9 is (0,2): walking up triples you.
• 6 is (1,1), 36 is (2,2), 12 is (2,1).

The grid isn’t a decoration. The position encodes the number, and — as we’re about to see — the operations you do to numbers become motions on the grid.

Multiplying is just adding the arrows

What happens to the coordinates when you multiply two numbers?

Multiplying means combining their prime factorizations — you pour both bags of primes into one bag. Powers of 2 add to powers of 2, powers of 3 add to powers of 3, each prime minding its own axis:

In coordinates that is the plainest thing in the world: add the vectors componentwise. Multiplication of numbers is addition of their coordinate vectors.

Take 12 = (2,1) and 18 = (1,2). Lay the second arrow tip-to-tail on the first, exactly the way you added vectors in school, and you land on (3,3) — which is 2³·3³ = 216. And indeed 12 × 18 = 216.

Everything you know about vector addition now says something about arithmetic.

Division is subtraction. To do 216 ÷ 18, walk the 216 arrow backwards by 18: (3,3) − (1,2) = (2,1) = 12. The catch: you only land on a whole number if no coordinate goes negative — which is exactly the condition that 18 divides 216. A negative coordinate would be a fraction.

Raising to a power is scaling. Squaring doubles every exponent, so it just stretches the arrow to twice its length: 12² = (4,2) = 144. Cubing scales it by 3, and so on.

And 1 is the origin — multiplying by 1 moves you nowhere, the zero vector.

The whole multiplicative world has quietly become a vector space over the exponents.

gcd and lcm are two corners of a box

This is the part that made me actually grin the first time I drew it.

Take two numbers and look at the smallest rectangle that has both of them as corners. Its other two corners — the bottom-left and the top-right — are the greatest common divisor and the least common multiple. For free.

Why? Go back to coordinates. A number dd divides n exactly when d sits “below and to the left” of n on every axis — every exponent of d is ≤ the matching exponent of n (you can’t have more 2s than n has and still divide it). So:

• the gcd of two numbers takes, on each axis, the smaller of the two exponents — the largest box that still fits underneath both. That’s the bottom-left corner.
• the lcm takes, on each axis, the larger exponent — the smallest box that covers both. That’s the top-right corner.

For 24 = (3,1) and 18 = (1,2):

• gcd: take the min on each axis → (1,1) = 2·3 = 6.
• lcm: take the max on each axis → (3,2) = 8·9 = 72.

Two famous facts fall out of the picture with no algebra at all.

First, gcd(a, b) × lcm(a, b) = a × b. On each axis, min + max is just the two original exponents added together (the smaller plus the larger is the same as the first plus the second). Adding coordinates is multiplying numbers — so the bottom-left corner times the top-right corner equals the two original numbers multiplied. Check it: 6 × 72 = 432 = 24 × 18. ✓

Second, two numbers are coprime exactly when their box is pinned to the axes — when the bottom-left corner is the origin (0,0) = 1. No shared axis, no common factor.

Adding a third prime is just adding a third axis

Nothing about this stops at two primes. A number whose factors are 2, 3, and 5 needs a third axis — the power of 5 — and now lives in 3-D space. The number 30 = 2·3·5 sits at (1,1,1), the unit corner pointing equally along all three primes:

Everything we just did still holds, verbatim. Multiply two numbers? Add their 3-D arrows. gcd and lcm? The min-corner and max-corner of the box they span, now a box in three dimensions. And it keeps going: a fourth prime is a fourth axis, a fifth prime a fifth axis, on and on. We just can’t draw past three — which is why the whole lattice is genuinely infinite-dimensional, with the primes as a basis. Every integer is one lattice point, reached by a unique combination of prime-axis steps.

So far we’ve stayed safely on the whole-number grid. But dividing and taking roots are about to carry us right off it — and each step we take outside rediscovers a famous kind of number.

Step one: negative coordinates are the fractions

We saw that dividing means walking an arrow backwards. So what happens when we divide too far — past the origin, into negative coordinates?

Take 6 = (1,1) and divide by 4 = (2,0). Subtracting walks us two steps left in the powers-of-2 direction, straight across the axis:

That negative coordinate isn’t a glitch. Allow negative exponents and the lattice fills in with every rational number.

Step two: off the grid entirely — the irrationals

Negative coordinates still land you on grid points. To leave the grid completely, we don’t subtract — we scale by a fraction.

This is the exact step the ancient Greeks refused to take. The Pythagoreans built their world on a single faith: that every quantity is a ratio of whole numbers — every number sitting, in our picture, on a tidy grid point. So when one of them (legend names Hippasus) proved that √2 is no such ratio — that it slips forever between the grid points and can never be pinned to one — it landed not as a theorem but as a scandal. The story even has him drowned at sea for letting it out. They did not want to leave the grid. We’re about to anyway :)

Remember that squaring doubled an arrow’s length. Running that backwards, a square root halves it. Take the square root of 6 = (1,1) and you scale the arrow by ½:

That point sits between the grid lines — at the exact midpoint of the diagonal out to 6, and there is no lattice crossing there. A point with non-integer coordinates is an irrational number.

√6 isn’t a failure of the system. It’s the system telling you, geometrically, why √6 can’t be a fraction: a fraction must land on a rational-coordinate point, and ½ along two prime axes at once never resolves to one.

And what about π, or e? They’re here too — they just sit at coordinates that are themselves irrational. To reach √6 we divided the exponent by a whole number (a square root), landing on the rational coordinate ½. But you can’t reach π by dividing by any whole or fractional number — you have to raise to an irrational power. Since π = 2^{1.6515…}, it sits on the 2-axis at height log₂π, and that height is provably irrational (if it were a fraction p/q, then 2ᵖ = πᵠ, impossible). So the hierarchy is beautifully clean:

• integer coordinates → whole numbers,
• negative coordinates included → fractions,
• rational coordinates → algebraic irrationals like √6, reached by roots,
• irrational coordinates → everything else, reached by irrational powers — where π, e, and the transcendentals live.

Nothing is missing from the space. The “wilder” a number is, the stranger the coordinates you need to pin it down. (A gem in passing: by the Gelfond–Schneider theorem, even an algebraic-but-irrational coordinate is enough to force transcendence — the point √2 on the 2-axis is 2^√2, a proven transcendental number.)

Step three: let the exponents be anything — and a number becomes a plane

Here’s where it gets strange in the best way. We’ve been treating 6 as the point (1,1) — the recipe 2 · 3. But is that the only way to reach the value 6?

What if we let the exponents be any real numbers? Then we could ask for pure powers:

Every one of these is a legitimate (a, b) whose value is exactly 6.

Collect them all, and they don’t form a point — they form a straight line:

The condition “2^a ⋅ 3^b = 6” is, after taking logs, just a⋅ln2 + b⋅ln3 = ln6 — the equation of a line. The integer recipe (1,1) is merely the one spot where that line happens to cross a grid point.

And with three primes the line becomes a plane — the set of all (a, b, c) with 2^a ⋅ 3^b ⋅ 5^c = 6, tilted through the 2-3-5 space and passing through the lone integer recipe at (1,1,0):

With all the primes, that plane becomes a hyperplane slicing through the infinite-dimensional space. That’s the final twist of the picture: once you free the exponents, a single number stops being a dot and reveals itself as an entire flat sheet of equivalent constructions — and the ordinary integer you started with is just where that sheet pierces the lattice.

The hidden ruler: how far apart are these sheets?

Here’s the prettiest part, and it’s exactly what makes the value a number rather than just a sheet. All of these “value = N” sheets are parallel — they share one normal direction, the log-prime vector

And there’s a clean identity behind it. For a point x=(a,b,c,…), the dot product with that normal is

So the sheet for value 1 passes through the origin (ln1 = 0), and as you push the sheet straight out along the normal, its value climbs: e^c where c is how far you’ve gone. Let it sweep outward through the 2–3–5 lattice and watch which integers it sweeps across — each one lights up at the exact moment the plane reaches its grid point:

And here’s the payoff that brings us all the way back to the start. The plane sails right past 7 and 11 without ever touching a point — because 7 and 11 are primes, and there is no point for them anywhere in the 2–3–5 space. They demand axes we didn’t draw. That’s the fundamental theorem of arithmetic staring back at us geometrically: a prime is a number with no grid point in any lower-prime space — it forces a brand-new dimension into being. Every other number up to 12 is just a corner reachable from 2, 3, and 5.

The gap is the same story a dimension down or up. The “value = N” sheet always sits at log-height lnN along the normal from the origin sheet, so 1, 6, 36, 216 — equally spaced multiplicatively — land as equally spaced parallel sheets:

Bigger numbers are farther out; multiplying numbers adds these heights (because logs add). The logarithm of a number is, quite literally, how far along the normal its sheet sits.

One honest footnote on the word “distance.” The exact, well-behaved quantity is the dot product n⋅x = lnN — a log-value, not a length. The true perpendicular distance is lnN divided by the length of the normal vector,

— the square root of the summed squares of the log-primes, one term per prime (not the count of primes). In the 2-D picture that length is √[(ln2)^2+(ln3)^2]​≈1.30, so the value-6 line sits about 1.38 from the origin, not 1.79. It’s a fixed constant in any finite picture, so the sheets stay equally spaced — just rescaled. But once you include all the primes, the sum under the root grows without bound, “Euclidean distance” collapses, and the log-value n⋅x is the only honest ruler.

One infinite sheet, one grid point

Now stop and really look at that value-6 sheet. It is infinite. It runs forever in every direction, slicing through an infinite-dimensional lattice that is itself packed with grid points everywhere you look — billions upon billions of them, out to infinity.

And it touches exactly one of them.

Not two. Not zero. Not a scattering of them every so often. The entire boundless sheet threads through the single point (1, 1, 0, …) — the factorization 2 · 3 — and misses every other grid point in the whole infinite universe of integers, forever.

That is unique factorization wearing its geometric face. A second grid point on this sheet would be a second whole-number recipe for 6 — a different bag of primes multiplying to the same value — and we know there isn’t one. So the sheet simply cannot graze another corner of the lattice. The uniqueness in the fundamental theorem is the statement that each value-sheet kisses the grid exactly once.

And now look back at the sheet sweeping outward. Every time it lit up an integer, it was kissing that number’s one grid point — so factoring a number is nothing more than knowing, in advance, which corner of the lattice its sheet is about to kiss. “Factor 60” means “predict the point the value-60 sheet will touch”: (2, 1, 1). And when the sheet glided past 7 and 11 touching nothing at all, that was the discovery of a prime — a sheet that kisses no point in the space it's allowed to reach.

And there’s something almost magical about why it can miss everything else so cleanly. The sheet is tilted by the log-primes ln2,ln3,ln5,…, and those numbers are incommensurable — no whole-number combination of them ever equals another (except the trivial, identical one). The tilt is “perfectly irrational” with respect to the grid. A sheet tilted at a rational angle would either miss the lattice entirely or, the instant it caught one point, catch an infinite row of them. This sheet does neither. It finds its one point and threads the infinite needle past all the rest — and it manages that trick precisely because of the primes underneath it.

That’s the picture I wanted to plant: the fundamental theorem of arithmetic isn’t only telling you that numbers factor. It’s handing you a coordinate system — and in it, integers are lattice points, fractions fill the quadrants, irrationals slip between the lines, and every value is a whole parallel sheet whose distance from the origin is its own logarithm.

The arithmetic you’ve done your whole life was geometry all along.

Imagine a flock of sheep scattered across a hillside, each one wearing a GPS collar. Every collar reports three numbers: how far east the sheep is, how far north, and its elevation. The sheep don’t move much — they just stand there, grazing.

From those readings alone, can you reconstruct the shape of the hill?

That’s the whole problem of regression in one picture. You have noisy measurements at scattered locations, and you want a model of the smooth surface underneath. The catch is that the tools most people reach for are surprisingly bad at it — and the tools that are good at it are black boxes you’re not allowed to use where we need to explain our decisions — like a bank or insurance company.

Let me show you what I mean, and then show you a model that finally gets to have both.

The usual suspects, and how they fail

Here’s every common regressor’s best attempt at recovering our hill. Each one was tuned to its best hyper-parameters on a held-out validation split, plus our new contestant.

Look at what happens:

A decision tree (R² = 0.84) turns a smooth hill into a ziggurat. Trees split the world into axis-aligned boxes and predict one flat value per box, so the surface comes out as a staircase. Two sheep standing one step apart on the hillside can get wildly different predicted elevations, just because a split happened to fall between them.

Linear regression (R² = 0.42) does something even worse: it tilts a flat plane through the data and calls it a day. A hill is not a plane. No amount of fitting fixes that — the model simply can’t represent a bump.

A random forest (R² = 0.94) and LightGBM (R² = 0.93) finally get the accuracy right. But notice the price: a random forest is 300 trees averaged together, and LightGBM is a few hundred boosted trees. There’s no single flowchart you can point at. And if you slice through their surfaces, you can see they’re not even smooth — they’re accurate but slightly jagged.

The fuzzy decision tree

A fuzzy decision tree keeps the thing that makes a tree understandable — a series of yes/no questions about your features — but throws out the part that makes it ugly. Instead of a hard split (“is elevation-east greater than 0.3? go left, else go right”), every split is a smooth transition. Near the boundary, a sample partially belongs to both sides, weighted by how far it sits from the threshold.

The consequence is exactly what you want on a hillside: neighbouring locations get neighbouring predictions. There are no cliffs between leaves. The prediction surface is continuous, and it bends instead of stepping.

And it’s still one tree, so you’d expect it to be weaker than a forest, but it isn’t:

One fuzzy tree beats the best random forest. On simpler problems it keeps pace with LightGBM too.

“But a single deep tree always overfits”

That’s the usual reason nobody trains one deep tree and trusts it. So let’s test it directly — grow the fuzzy tree deeper and deeper, and watch it against a random forest and LightGBM at the same depth:

Two things stand out. First, at matched depth the single fuzzy tree (blue) sits right on top of the 300-tree forest’s accuracy — it just gets there with one tree instead of three hundred. Second, and more importantly, it doesn’t blow up as it gets deeper. The fuzziness is the regularizer. Where a normal decision tree would carve the training data into one-sheep boxes and memorize the noise, the fuzzy memberships smear each sample across neighbouring leaves, so depth buys you expressiveness without buying you overfitting. The train-to-test gap stays small even with a hundred-plus leaves.

In code, that’s the whole model:

pip install FuzzyDTree

from fuzzydtree import FuzzyTreeRegressor

model = FuzzyTreeRegressor(
    max_depth=12, min_samples_leaf=1, 
    include_hard_splits=False, margin_min_scale=0.4
)
model.fit(X_train, y_train)

model.predict(X_test)            # smooth predictions
model.score(X_test, y_test)      # R²

Two of those arguments are worth understanding, because they decide whether the tree comes out smooth or staircased:

• include_hard_splits=False removes the option of a hard, zero-width step at a split. It’s important because a hard split always separates the training data most sharply, so if you leave it available the tree tends to pick it everywhere and you're back to a staircase-ish estimate. Turning it off forces every boundary to be a genuine ramp — which is the whole point on a smooth hill. Not every problem requires it, and the comprehensive test shows that keeping it True as default helps.
• margin_min_scale=0.4 sets how soft those ramps must be: the narrowest allowed transition is 0.4× the feature's local spread. Raise it for smoother surfaces, lower it toward the default when some features genuinely have sharp thresholds.

The flip side: on data with crisp, step-like signals — think “has this borrower ever been 90+ days late?” — you want hard splits, so you leave include_hard_splits=True (the default). The knobs let the tree match the shape of your problem; the README has the full guidance.

It’s a scikit-learn estimator — fit, predict, score, get_params, the works. Drop it in wherever you'd put a DecisionTreeRegressor.

Why a bank should care

Now replace “elevation on a hill” with “probability this customer defaults,” and the sheep with loan applicants. This is where the smoothness-and-explainability combo stops being an aesthetic preference and becomes a business and regulatory requirement.

Banks are, in practice, forced to use simple models — logistic regression or a shallow decision tree — for any decision that affects a customer, because they have to explain it. “Your loan was declined” has to come with reasons, and those reasons have to be honest reflections of what the model actually did, not a story bolted on afterward. That rules out random forests and gradient boosting for a lot of high-stakes decisions, and it leaves a painful amount of accuracy on the table.

A fuzzy decision tree breaks that trade-off. It’s a single tree, so every prediction is just a weighted walk from the root to the leaves — and you can read that walk off as a waterfall of feature contributions, built in. No SHAP, no surrogate model, no post-hoc approximation. The explanation is the model.

import matplotlib.pyplot as plt
from fuzzydtree import FuzzyTreeClassifier

clf = FuzzyTreeClassifier(max_depth=5)
clf.fit(X_applicants, defaulted)

# Why did THIS applicant get THIS decision?
fig, ax = plt.subplots(figsize=(9, 5))
clf.plot_log_odds_contributions(applicant, feature_names=feature_names, ax=ax)
plt.tight_layout()
plt.show()

On the left, an applicant the model declines: a high debt-to-income ratio and high revolving utilization push the default odds up, and not even a decent income pulls them back enough. On the right, an approved applicant: strong income dominates, and nothing else moves the needle much. (This example uses a synthetic credit dataset for illustration, but the mechanics are identical on real ones — e.g. Kaggle’s “Give Me Some Credit” or the German Credit data.)

A loan officer can read those two charts. A regulator can read them. The customer can read them. And the model behind them is as capable as the random forest the bank wasn’t allowed to deploy.

How it compares on ordinary benchmarks

The honest benchmark story is not “one fuzzy tree always beats a forest.” It usually doesn’t. If your only goal is the highest possible predictive score, a random forest or a gradient-boosted ensemble is still often the better default: they get to average or boost many trees, and that extra capacity matters.

The more interesting comparison is against a single hard decision tree, because that is the model class people choose when they need something compact, inspectable, and easy to explain. On that comparison, the fuzzy tree is often a large upgrade.

In the regression benchmarks, random forests are frequently ahead of one fuzzy tree, and LightGBM is often ahead of both. But the fuzzy tree consistently behaves like a much stronger single-tree model than DecisionTreeRegressor. The hard tree tends to plateau early or degrade as depth increases; the fuzzy tree keeps improving on several datasets because soft membership gives it a smoother bias without turning it into an ensemble.

The classification picture is similar. A random forest or LightGBM is usually the accuracy-first choice, especially when the data has many interactions or noisy boundaries. But the fuzzy classifier can easily beat a lone hard decision tree while remaining a single model whose predictions can be traced and explained directly.

So the practical positioning is:

• If you are free to deploy a black-box ensemble and only care about score, use the forest or booster.
• If you need a single transparent model, a fuzzy tree is often a better baseline than a plain decision tree.
• If you need smooth predictions or stable explanations near split boundaries, the fuzzy tree gives you something ordinary trees simply cannot represent.

The point

For years the implicit deal in applied machine learning has been: if you want accuracy, accept a black box; if you want to explain yourself, accept a worse model. Linear models and lone decision trees on one side; random forests and gradient boosting on the other.

The sheep on the hill show that deal isn’t necessary. A single fuzzy decision tree:

• recovers a smooth surface where a plain tree gives a staircase and a linear model gives a plane,
• can approach forest-level accuracy on smooth problems — with one tree,
• often beats a plain single decision tree on tabular benchmarks,
• doesn’t overfit as it deepens, because the fuzziness regularizes it, and
• explains every prediction with a built-in waterfall, no SHAP required.

You’re no longer forced to pick a black box just to get good predictions — for regression or classification.

pip install FuzzyDTree

• Docs: https://guyko81.github.io/fuzzydtree/
• GitHub: https://github.com/guyko81/fuzzydtree
• PyPI: https://pypi.org/project/fuzzydtree/

The short version: I believe Mythos combines a looped Transformer — the same block of layers running multiple times, giving the model deeper reasoning without more parameters — with an Engram-style memory layer that handles fact retrieval in constant time, freeing every loop iteration for pure reasoning. The rest of this article is the evidence trail.

On March 26, 2026, Fortune broke the story: a draft blog post in an unsecured Anthropic data store revealed a model called Claude Mythos — described internally as a “step change” in capabilities. Two weeks later, Anthropic officially announced Mythos Preview through Project Glasswing, a cybersecurity alliance with Amazon, Apple, Google, Microsoft, Nvidia, CrowdStrike, and others.

Leaked internal documents call it the “Capybara” tier — a new model class above Opus, “larger and more intelligent” than anything before it.

The benchmarks show why. From the system card (Table 6.3.A):

Look at GraphWalks. Mythos more than doubles Opus 4.6 and nearly quadruples GPT-5.4 on a benchmark that requires coherent reasoning across up to a million tokens of graph traversal.

And in cybersecurity, it autonomously found thousands of zero-day vulnerabilities — including a 27-year-old bug in OpenBSD and a 17-year-old remote code execution flaw in FreeBSD — bugs that survived decades of human review and millions of automated fuzzing tests.

Something fundamental changed. But raw scale alone doesn’t explain it — GPT-5.4 is massive too, yet Mythos outperforms it by 4x on long-context graph reasoning. The question: what architectural innovation makes Mythos different?

The answer, I believe, lies in a quiet revolution unfolding in the research literature since late 2025.

The Looped Transformer Revolution

Starting in 2025, a series of papers converged on an alternative: what if you ran the same layer block multiple times, like a loop? The model doesn’t need more parameters — it needs more iterations. Each pass refines the internal representation in high-dimensional latent space, and the model only projects back to tokens at the very end.

The idea has an ancestor.

The Universal Transformer (2018) — The Seed

Google’s Universal Transformer (Dehghani et al., ICLR 2019) was the first to bring recursion into Transformers. A single block repeated T times, with an Adaptive Computation Time (ACT) mechanism that learns when to stop:

Same function, every step. A fixed-point iteration that refines until it converges.

It was ahead of its time — the scaling laws of the era rewarded more parameters, not more iterations. But the idea never died.

HRM (June 2025) — The Viral Moment

The Hierarchical Reasoning Model (HRM) changed everything. Two recurrent modules — one for slow strategic planning, one for fast detailed computation — inspired by the prefrontal cortex vs. lower brain regions. With just 27M parameters and 1,000 training samples, it scored 40.3% on ARC-AGI-1.

It went viral. Over 4 million views on social media. The ARC Prize Foundation’s analysis later showed much of the gain came from the outer loop refinement process rather than the hierarchy itself — but the genie was out of the bottle. The community was paying attention to recursive architectures again.

TRM (October 2025) — The Simplification

Samsung SAIL’s Tiny Recursive Model stripped HRM down to its essence: one module, two roles. A 2-layer, 7M-parameter network alternating between “thinking” (updating the latent scratchpad with the input) and “answering” (updating the solution without the input). It hit 45% on ARC-AGI-1 — outperforming models 10,000x its size — with full backpropagation through all recursive steps.

The message was clear: recursive depth beats parameter count for reasoning.

ByteDance Enters — And Solves the Scaling Problem

The small recursive models proved the concept, but they had a fatal flaw for production use: loops run sequentially. If you have L loops, every token takes L times longer to generate. Fine for a 7M research model. Catastrophic for a frontier-scale production system.

ByteDance published two papers in October 2025 that solved this.

Parallel Loop Transformer / PLT — Zero Extra Latency

The Parallel Loop Transformer introduced Cross-Loop Parallelism (CLP): instead of waiting for all loops of one token to finish, you compute different loops of different tokens simultaneously:

• Token i, Loop 1
• Token i-1, Loop 2
• Token i-2, Loop 3
• …all in the same decoding step.

This works because LLM decoding is memory-bound, not compute-bound — most of the GPU is sitting idle anyway. The parallel FLOPs are essentially free.

Combined with Gated Sliding-Window Attention (sharing the first loop’s KV cache globally), PLT achieves the accuracy of a traditional looped model with virtually zero extra latency.

Sound like exactly what you’d need to deploy a looped model at frontier scale?

Ouro / LoopLM — The Large-Scale Proof

Ouro (“Ouroboros”) is the first large-scale, pre-trained looped language model. A standard decoder-only Transformer (RoPE, SwiGLU, sandwich norm) whose entire block runs T=4 times. The latent state never collapses into token space between loops — the model thinks in high-dimensional latent space and only projects to tokens at the very end.

The headline result: Ouro-1.4B with 4 recurrent steps matches 4B-parameter standard transformers.

But here’s the result that matters most for our theory. The ByteDance team ran controlled experiments and found:

The advantage comes from superior knowledge manipulation — especially on tasks requiring multi-step reasoning and graph search.

Now look at that Mythos benchmark table again. GraphWalks BFS: 80%. Opus 4.6: 38.7%. GPT-5.4: 21.4%.

If the Ouro findings generalize — that looping specifically improves knowledge manipulation over long, structured contexts — then a looped architecture would predict exactly this pattern: a disproportionate jump on graph traversal, the benchmark most dependent on iterative multi-hop reasoning.

The Missing Piece: Where Does the Knowledge Live?

Looped architectures give you deeper reasoning, and PLT makes it deployable. But Ouro’s blind spot remains: ~2 bits per parameter, no matter how many loops. Loops help you use what you know, but they don’t help you know more.

The problem runs deeper. A conventional Transformer stores everything in its weights — facts, entities, idioms, all tangled together with the computation that’s supposed to reason over them. When the model encounters “Diana, Princess of Wales,” it doesn’t just look up who that is. It burns computation cycles in the early layers assembling those three tokens into a coherent entity — rebuilding a lookup table from scratch, every single forward pass. That’s a chess grandmaster re-deriving the rules of chess before every move. The knowledge is there, but it’s encoded in the wrong format — mixed into the same weights that do the thinking.

For a model like Mythos that needs to reason over vast codebases, chain obscure CVEs together, and recall the specifics of a 27-year-old OpenBSD system call — you need a different approach entirely.

DeepSeek Engram (January 2026) — The Memory Primitive

DeepSeek’s Engram (code) proposes exactly that: pull the static knowledge out into a dedicated O(1) lookup module, and let the Transformer layers focus on reasoning. Deterministic hash functions map N-grams to embedding tables — no attention, no FFN, just a constant-time table lookup. A learned gating mechanism decides when to trust the result and when to ignore it, since hash collisions mean the lookup sometimes returns garbage.

The key infrastructure trick: because the lookups are deterministic and read-only during inference, the entire memory module can live on CPU RAM with negligible latency. The model’s factual capacity is no longer bounded by GPU memory — you can keep scaling the lookup tables on cheap system memory while the expensive GPU stays dedicated to active computation.

In the Engram-27B model, DeepSeek traded 17 MoE experts for a 5.7 billion parameter embedding module — and performance went up across the board. The lesson: when you give a model a proper place to store facts, the layers you free up become better at everything else.

For a looped model, this matters even more. Every loop iteration that used to be partly wasted on recalling facts can now focus entirely on what loops are good at: connecting dots, chaining inferences, exploring paths. One provides the facts. The other provides the thinking. For a model that chains KASLR bypasses, heap sprays, and ROP gadgets into working exploits in under half a day, this pairing starts to look less like speculation and more like a requirement.

Putting It Together: The Mythos Theory

Anthropic has said nothing about what’s under the hood. This is a theory — but the circumstantial case is hard to ignore.

ByteDance showed you can loop for depth and parallelize the loops. DeepSeek showed you can pull knowledge into a separate memory module. Maybe Anthropic put the pieces together.

Mythos scores higher in exactly the places you’d expect looping to matter: 4x on graph traversal over GPT-5.4, disproportionate jumps in multi-step exploit chaining and long-context coherence. Its cybersecurity capabilities “emerged as a downstream consequence” of general improvements — not domain-specific training.

But Mythos also knows more. MMMLU up to 92.7% from Opus 4.6’s 91.1%, GPQA Diamond to 94.5% from 91.3%. Ouro showed looping can’t push knowledge storage past ~2 bits per parameter — so where are the extra facts coming from? That’s the Engram-shaped hole.

Maybe Anthropic found a completely different path. But they gave Mythos a new tier name because the old ones didn’t fit. I think it’s not a larger Transformer, a bigger Opus. I think it’s a new era.

• Universal Transformers — Dehghani et al., 2018
• Hierarchical Reasoning Model — Wang et al., 2025
• Less is More: Recursive Reasoning with Tiny Networks (TRM) — Jolicoeur-Martineau et al., 2025
• Parallel Loop Transformer — ByteDance, 2025
• Scaling Latent Reasoning via Looped Language Models (Ouro) — ByteDance, 2025
• Conditional Memory via Scalable Lookup (Engram) — DeepSeek, 2026
• Project Glasswing — Anthropic, 2026
• Claude Mythos Preview System Card — Anthropic, 2026
• Fortune: Anthropic ‘Mythos’ AI model representing ‘step change’

When your model predicts “tomorrow’s sales will be 4,200 units,” what does that actually mean? Is it confident? Could it be 3,000 or 6,000 with equal likelihood? Is the distribution symmetric, or is there a long tail of risk on one side?

A point prediction tells you none of this. And yet, most regression models in production output exactly one number per sample.

This matters: if you’re setting inventory levels, pricing insurance, or allocating resources, the shape of the uncertainty around a prediction is often more valuable than the prediction itself. A bimodal distribution (two likely outcomes) requires a fundamentally different decision than a tight unimodal one, even if the mean is the same.

DistributionRegressor is a Python package that predicts full probability distributions p(y|x) using a single LightGBM model. It’s nonparametric, fast, and requires no distributional assumptions. In this article, I’ll explain how it works and why it takes a different approach than existing methods.

pip install distribution-regressor

The problem with parametric distributional regression

The most well-known approach to distributional regression is NGBoost (Duan et al., 2020), which uses natural gradient boosting to fit parameters of a chosen distribution family — Normal, LogNormal, etc.

NGBoost is elegant, but it has a fundamental constraint: you must choose the parametric family in advance. If you choose Normal but the true conditional distribution is bimodal, every prediction will be wrong in a way that no amount of training data can fix. The model literally cannot represent the true distribution.

This isn’t a theoretical concern. Real-world conditional distributions are often:

• Multimodal — commute times (traffic vs. no traffic), house values (neighborhood mix)
• Zero-inflated — retail sales (store open vs. closed), insurance claims
• Heavily skewed with complex tails — medical costs, financial returns

You could try fitting a mixture model, but then you need to choose the number of components. You could use a more flexible parametric family, but you’re always one assumption away from misspecification.

A different approach: learn the CDF directly

DistributionRegressor sidesteps parametric assumptions entirely. Instead of fitting distribution parameters, it learns the conditional cumulative distribution function (CDF) directly:

The key insight is that this is a binary classification problem. For any threshold t and any training sample (x_i, y_i), we know exactly whether y_i falls below τ:

This is just a binary label — and LightGBM is very good at binary classification.

The algorithm in four steps

1. Create a grid. Discretize the target range into n_bins threshold points τ_1,τ_2,…,τ_K.

2. Create binary targets. For each training sample and each threshold, compute z_ij​=1{y_i​≤τ_j​}.

3. Expand the dataset. Cross every sample with every threshold point. A dataset of N samples becomes N×K rows, each with the original features plus the threshold value as an additional feature.

4. Train a single LightGBM model with cross-entropy loss and a monotone increasing constraint on the threshold feature. This constraint is critical — it ensures the model outputs a valid CDF (non-decreasing in τ).

At inference, the model predicts F(τ∣x) for all grid points, and the probability mass function (PMF) is recovered by differencing the CDF.

That’s it. One model. No custom loss functions, no natural gradients, no distribution parameters. The monotonicity constraint does the heavy lifting — it forces the output to be a valid CDF regardless of the input.

Why this works surprisingly well

The CDF formulation has several elegant properties:

Every training sample provides dense signal. A single sample (x_i​,y_i​) with K grid points generates K training examples: all thresholds below yi_​ get target 0, all thresholds above get target 1. There’s no ambiguity — the labels are exact.

Monotonicity is baked in. LightGBM’s native monotone constraints ensure the predicted CDF is non-decreasing in the threshold, so you always get a valid probability distribution.

Logistic loss gives calibrated probabilities. Cross-entropy loss naturally encourages the model to output well-calibrated probabilities, which is exactly what we want for CDF values.

Arbitrary distribution shapes are free. Multimodal? Skewed? Zero-inflated? The model doesn’t care — it just learns where the CDF steps are.

Show me the code

from distribution_regressor import DistributionRegressor

# Train - that's it
model = DistributionRegressor(
    n_bins=100,
    n_estimators=1000,
    learning_rate=0.01,
)
model.fit(X_train, y_train)

# Point predictions
y_mean = model.predict(X_test)          # E[Y|X] (mean)
y_mode = model.predict_mode(X_test)     # peak of density

# Uncertainty quantification
intervals = model.predict_interval(X_test, confidence=0.9)
std = model.predict_std(X_test)

# Full distribution
grids, dists, _ = model.predict_distribution(X_test)

# Proper scoring
nll = model.nll(X_test, y_test)  # negative log-likelihood

The API follows scikit-learn conventions. It also provides pdf(), cdf(), ppf() (inverse CDF), and logpdf() for density evaluation at arbitrary points.

Real-world example: retail sales forecasting

Let’s look at a concrete example using the Rossmann Store Sales dataset. We train a model to predict daily sales for a single store using calendar features, promotions, and holiday indicators.

model = DistributionRegressor(
    n_bins=100,
    n_estimators=1000,
    learning_rate=0.01,
)
model.fit(train[features], train['Sales'])

Time series with uncertainty bands

The model produces prediction intervals that naturally adapt to context. Promotional days show wider uncertainty bands. Sundays (when the store is closed) collapse to zero. The model learns all of this from the data — no special handling required.

Individual predicted densities

Each test sample gets a full predicted density. The closed-day prediction is a sharp spike at zero. Promotional days show wider spreads. Regular weekdays are tighter. The red dashed line is the actual sales value.

These are not Gaussian distributions forced onto the data. They’re whatever shape the data requires.

Comparison with NGBoost on California Housing

To benchmark against NGBoost, we use the California Housing dataset from scikit-learn (20,640 samples, 8 features). House value distributions are often multimodal and bounded — properties that violate simple parametric assumptions.

DistributionRegressor outperforms NGBoost by 0.42 nats — without any distributional assumptions. NGBoost’s Normal distribution forces every prediction into a symmetric bell curve, which is a poor fit for bounded, multimodal house value distributions. (Target is median house value in $100k, sklearn default scale.)

DistributionRegressor (blue) captures complex conditional densities that NGBoost’s Normal (green dashed) simply cannot represent. The parametric model is forced into a symmetric shape regardless of the true data distribution.

The nonparametric approach wins on both aggregate NLL and individual prediction quality. If your downstream decision depends on the shape of the distribution (not just its mean or variance), the nonparametric approach gives you information that the parametric one discards.

What about memory?

The main practical consideration is memory. The training set expansion — N samples times K grid points — means a 100k-sample dataset with 100 bins becomes 10M rows. This is manageable for most problems, but can be a constraint for very large datasets.

For those cases, there’s a Monte Carlo training mode (monte_carlo_training=True) that samples a subset of grid points per training sample, reducing peak memory by 5-20x. The mc_resample_freq parameter resamples fresh grid points every N trees during training, so different groups of trees see different thresholds — analogous to mini-batch SGD over the grid dimension. In our benchmarks, MC with resampling produced distributions close to the full grid at a fraction of the memory cost.

Key design decisions

A few choices that make this work in practice:

Single model, not an ensemble of binary classifiers. Earlier approaches (like the histogram-based methods) train a separate classifier for each threshold. DistributionRegressor uses the threshold as a feature in a single model, which lets it share information across thresholds and requires far fewer parameters.

Cross-entropy loss, not MSE. The binary targets z_ij​∈{0,1}, and we want calibrated probabilities. Cross-entropy is the natural choice and gives much better-calibrated CDFs than MSE.

Monotone constraints, not post-hoc correction. Rather than predicting unconstrained values and fixing monotonicity violations after the fact, the constraint is built into the LightGBM training. This produces smoother, more stable CDFs.

Output smoothing as a post-processing step. An optional Gaussian smoothing on the predicted CDF can be applied (and tuned) after training, without refitting. This is useful for smoothing stepped densities at lower grid resolutions, and can preserve discrete mass points (like exact zeros) through the atom_values parameter.

When to use this

DistributionRegressor is a good fit when:

• You need full conditional distributions, not just point predictions
• The true distribution shape is unknown or complex (multimodal, skewed, zero-inflated)
• You want a nonparametric approach with no distributional assumptions
• Your data is tabular (it’s LightGBM under the hood)
• You want a scikit-learn compatible API that drops into existing workflows

It’s not the right tool if you need distributions over structured outputs (sequences, images), if your data is very high-dimensional (use a neural density estimator), or if you have strong prior knowledge about the parametric family (NGBoost will be more accurate in that case).

Try it

pip install distribution-regressor

• Documentation: guyko81.github.io/DistributionRegressor
• GitHub: github.com/guyko81/DistributionRegressor
• PyPI: pypi.org/project/distribution-regressor

The package is MIT licensed and the API is stable at v2.0. Contributions and feedback are welcome.

Gabor Gulyas