Two Hours of Experimental Mathematics

Stephen Wolfram leading a live experiment at MoMath

A Talk, a Performance… a Live Experiment

“In the next hour I’m going to try to make a new discovery in mathematics.” So I began a few days ago at two different hour-long Math Encounters events at the National Museum of Mathematics (“MoMath”) in New York City. I’ve been a trustee of the museum since before it opened in 2012, and I was looking forward to spending a couple of hours trying to “make some math” there with a couple of eclectic audiences from kids to retirees.

People usually assume that new discoveries aren’t things one can ever see being made in real time. But the wonderful thing about the computational tools I’ve spent decades building is that they make it so fast to implement ideas that it becomes realistic to make discoveries as a kind of real-time performance art.
Try the experiments for yourself in the Wolfram Open Cloud »

But mathematics is an old field. Haven’t all the “easy” discoveries already been made? Absolutely not! Mathematics has progressed along definite lines, steadily adding theorems about all sorts of things. Many great mathematicians (Gauss, Ramanujan, etc.) have done experimental mathematics to find out what’s true. But in general, experimental mathematics hasn’t been pursued nearly as much as it could or should have been. And that means that there’s a huge amount of “low-hanging fruit” still to be picked—even if one’s only got a couple of hours to spend doing it.

Experiment #1

My rule for live experiments is that to keep everything fresh I think of the topic only a few minutes before I start. But since this was my first-ever such event at the museum, I thought I should have a topic that’s somehow a big one for me. So I decided it should be something related to cellular automata—which were the very first examples I explored in the multi-decade journey that led to A New Kind of Science.

While their setup is nice and easy to understand, cellular automata are fundamentally systems from the computational universe, not “mathematical” systems. But what I thought I’d do for my first experiment was to look at some cellular automata that somehow have a traditional mathematical interpretation.

After introducing cellular automata (and the Wolfram Language), I started off talking about Pascal’s triangle—formed by making each number to be the sum of left and right neighbors at each step. Here’s the code I wrote to make Pascal’s triangle (yes, replacing 0 by “” is a bit hacky, but it makes everything much easier to read):

NestList[RotateLeft[#] + RotateRight[#] &, CenterArray[{1}, 21, 0, 10]/. 0 -> "" // Grid

If one does the same thing mod 2, one gets a rather clear pattern:

NestList[Mod[RotateLeft[#] + RotateRight[#],2] &, CenterArray[{1}, 21, 0], 10]/. 0 -> "" // Grid

And one can think of this as a cellular automaton, with this rule:

RulePlot[CellularAutomaton[90]]

Here’s what happens if one runs this, starting from a single 1:

ArrayPlot[CellularAutomaton[90, {{1}, 0}, 50]]

And here’s the same result, from the “mathematical” code:

ArrayPlot[  NestList[Mod[RotateLeft[#] + RotateRight[#], 2] &,    CenterArray[{1}, 101, 0], 50]]

OK, so what happens if we change the math a bit? Instead of using mod 2, let’s use mod 5:

ArrayPlot[  NestList[Mod[RotateLeft[#] + RotateRight[#], 5] &,    CenterArray[{1}, 101, 0], 50]]

It’s still a regular pattern. But here was my idea for the experiment: explore what happens if the rule involves mathematical operations other than pure addition.

So what about multiplication? I was mindful of the fact that all the 0s in the initial conditions would tend to make a lot of 0s. So I thought: let’s try adding constants before doing the multiplication. And here’s the first thing I tried:

ArrayPlot[  NestList[Mod[(1 + RotateLeft[#])*(2 + RotateRight[#]), 5] &,    CenterArray[{1}, 101, 0], 50]]

I was pretty surprised. I wasn’t expecting anything that complicated. But, OK, I thought, let’s back off and try an even simpler rule: let’s use mod 3 instead of mod 5. (Mod 2 would already have been covered by my exhaustive study of the “elementary cellular automata”.)

Here’s the result I got:

ArrayPlot[  NestList[Mod[(1 + RotateLeft[#])*(2 + RotateRight[#]), 3] &,    CenterArray[{1}, 101, 0], 50]]

I immediately said, “I wonder how fast that pattern grows.” I guessed it might be a logarithm or a square root.

But before going on, I wanted to scope out what else was there in the space of rules like this. Just to check, I ran the mod 2 case. As expected, nothing interesting.

Table[ArrayPlot[   NestList[Mod[(a + RotateLeft[#])*(b + RotateRight[#]), 2] &,     CenterArray[{1}, 101, 0], 50]], {a, 0, 1}, {b, 0, 1}]

OK, now the mod 3 case:

Table[ArrayPlot[   NestList[Mod[(a + RotateLeft[#])*(b + RotateRight[#]), 3] &,     CenterArray[{1}, 101, 0], 50]], {a, 0, 2}, {b, 0, 2}]

An interesting little collection. But then it was time to analyze the growth of those patterns.

The first step, as suggested by someone in the audience, was just to rotate the list every step, to make the straight edge be vertical:

ArrayPlot[  NestList[RotateRight[     Mod[(1 + RotateLeft[#])*(2 + RotateRight[#]), 3]] &,    CenterArray[{1}, 100, 0], 50]]

Then we picked every other step, to get rid of the horizontal stripes:

ArrayPlot[  Take[NestList[    RotateRight[Mod[(1 + RotateLeft[#])*(2 + RotateRight[#]), 3]] &,     CenterArray[{1}, 200, 0], 200], 1 ;; -1 ;; 2]]

And—when in doubt—just run it longer, here for 3000 steps. Well, my guess about square root or logarithm was wrong: this looks roughly linear, albeit irregular.

ArrayPlot[  Take[NestList[    RotateRight[Mod[(1 + RotateLeft[#])*(2 + RotateRight[#]), 3]] &,     CenterArray[{1}, 3000, 0], 3000], 1 ;; -1 ;; 2]]

I was disappointed that this was so gray and hard to read. Trying colors didn’t help, though; the pattern is just rather sparse.

Well, then I tried to just plot the position of the right-hand edge. Here’s the code I came up with:

data = (First[#] - #) &[    Flatten[FirstPosition[Reverse[#], 1 | 2] & /@       Take[NestList[        RotateRight[          Mod[(1 + RotateLeft[#])*(2 + RotateRight[#]), 3]] &,         CenterArray[{1}, 2000, 0], 2000], 1 ;; -1 ;; 2]]];

Here’s a fit:

Fit[data, {1, x}, x]

OK, how can one get some better analysis? First, I took differences to see the growth at each step: always either 0 or 2 cells. Then I looked for runs of growth or no growth. And then I looked specifically for runs of growth, and saw how long the successive runs were.

runs = Length /@ Take[Split[Differences[data]], 1 ;; -1 ;; 2];

What is this? Being New York, there were lots of finance people in the audience—including in the front row a world expert on power laws. So the obvious question was, did the spikes have a power-law distribution of sizes? The results, based on the data I had, were inconclusive:

Histogram[runs]

But instead of looking further at this particular rule, I decided to take a quick look at the case of higher moduli. These are the results I got for mod 4:

Table[ArrayPlot[   NestList[Mod[(a + RotateLeft[#])*(b + RotateRight[#]), 4] &,     CenterArray[{1}, 100, 0], 50], PlotLabel -> {a, b}], {a, 0, 3}, {b,    0, 3}]

There was one that looked interesting here:

ArrayPlot[  NestList[Mod[(1 + RotateLeft[#])*(1 + RotateRight[#]), 4] &,    CenterArray[{1, 2, 1}, 100, 0], 50]]

Would it end up having lots of different possible structures? Trying it with random initial conditions made it look like it was never going to have anything other than repetitive behavior:

ArrayPlot[  NestList[Mod[(1 + RotateLeft[#])*(1 + RotateRight[#]), 4] &,    RandomInteger[2, 100], 50]]

Well, by this point our time was basically up. But it was hard to stop. I quickly tried the case of mod 5—and discovered all sorts of interesting behavior:

Table[ArrayPlot[   NestList[Mod[(a + RotateLeft[#])*(b + RotateRight[#]), 5] &,     CenterArray[{1, 2}, 100, 0], 50], PlotLabel -> {a, b}], {a, 0,    4}, {b, 0, 4}]

I just had to check out a couple of these. One that has an overall nested pattern, but with lots of complicated stuff going on “in the background”:

ArrayPlot[  NestList[Mod[(3 + RotateLeft[#])*(3 + RotateRight[#]), 5] &,    CenterArray[{1, 2}, 400, 0], 200]]

And one with a mixture of regular and irregular growth:

ArrayPlot[  Take[NestList[Mod[(2 + RotateLeft[#])*(4 + RotateRight[#]), 5] &,     CenterArray[{1, 2}, 2000, 0], 1000], 1 ;; -1 ;; 2]]

It was time to stop. But I was pretty satisfied. Live experiments are always risky. And we might have found nothing interesting. But instead we found something really interesting: rich and complex behavior based on iterating rules given by simple algebraic formulas. In a sense what we found is an example of a bridge between traditional mathematical constructs (like algebraic formulas), and pure computational systems, with arbitrary computational rules. In an hour we certainly didn’t finish—but we found a seed for all sorts of future research—on what we might call “MoMath Cellular Automata.”

Experiment #2

After a break, it was time for experiment #2. This time I decided to do something more related to numbers. I started by talking about reversal-addition systems—where at each step one adds a number to the number obtained by reversing its digits. I showed the result for base 10, starting from the number 123:

PadLeft[IntegerDigits[    NestList[FromDigits[Reverse[IntegerDigits[#]]] + # &, 123,      100]]] // ArrayPlot

Then I said, “Instead of reversing the digits, let’s just rotate them to the left. And let’s make the system simpler, by using base 2 instead of base 10.”

This was the sequence of numbers obtained, starting from 1:

NestList[FromDigits[RotateLeft[IntegerDigits[#, 2]], 2] + # &, 1, 10]

Someone asked whether it was a recognizable sequence. FindSequenceFunction didn’t think so:

FindSequenceFunction[%]

Then the question was, what’s the overall pattern? Here’s the result for 100 steps:

PadLeft[IntegerDigits[    NestList[FromDigits[RotateLeft[IntegerDigits[#, 2]], 2] + # &, 1,      100], 2]] // ArrayPlot" src="https://content.wolfram.com/sites/43/2017/03/21.png

It looks remarkably complex. And doing 1000 steps doesn’t make it look any simpler:

PadLeft[IntegerDigits[    NestList[FromDigits[RotateLeft[IntegerDigits[#, 2]], 2] + # &, 1,      1000], 2]] // ArrayPlot" src="https://content.wolfram.com/sites/43/2017/03/22.png

What about starting with something other than 1?

Table[PadLeft[    IntegerDigits[     NestList[FromDigits[RotateLeft[IntegerDigits[#, 2]], 2] + # &, n,       100], 2]] // ArrayPlot, {n, 10}]

All pretty similar. I wondered if rotating right, rather than left, would make a difference. It really didn’t:

Table[PadLeft[IntegerDigits[     NestList[FromDigits[RotateLeft[IntegerDigits[#, 2]], 2] + # &, n,       100], 2]] // ArrayPlot, {n,10}]

I thought maybe it’d be interesting to have a fixed number of digits, so I tried reducing mod 220, to keep only the last 20 digits:

Table[PadLeft[    IntegerDigits[     NestList[FromDigits[RotateRight[IntegerDigits[#, 2]], 2], 2^20] + # &, n,       100], 2]] // ArrayPlot, {n, 10}]

Table[FindTransientRepeat[NestList[Mod[FromDigits[RotateRight[IntegerDigits[#,2]],2+#,2^n]&,1,1000],4,{n,8}] // Column" src="https://content.wolfram.com/sites/43/2017/03/26.png

Then I decided to make complete transition graphs for all 2n states in each case. Curious-looking pictures, but not immediately illuminating.

Table[Labeled[   Graph[# ->        Mod[FromDigits[RotateRight[IntegerDigits[#, 2]], 2] + #,         2^n] & /@ Range[2^n]], n], {n, 2, 9}]

By now I was wondering: “Is there a still simpler system involving digit rotation that does something interesting?” I wondered what would happen if instead of adding in the original number at each step, I just multiplied by 2, and added some constant. This didn’t immediately lead to anything interesting:

Table[PadLeft[    IntegerDigits[     NestList[2 FromDigits[RotateLeft[IntegerDigits[#, 2]], 2] + n &,       1, 100], 2]] // ArrayPlot, {n, 6}]

So then I wondered about multiplying by 3:

Table[PadLeft[    IntegerDigits[     NestList[3 FromDigits[RotateLeft[IntegerDigits[#, 2]], 2] + n &,       1, 100], 2]] // ArrayPlot, {n, 7}]

Again, nothing too exciting. But—just to be complete—I thought I’d better run the experiment of looking at a sequence of other multipliers.

Table[Labeled[   PadLeft[IntegerDigits[      NestList[a FromDigits[RotateLeft[IntegerDigits[#, 2]], 2] + 1 &,        1, 100], 2]] // ArrayPlot, a], {a, 20}]

Similar behavior until—aha!—something weird and complicated happens when one gets to multiplier 13.

There was an immediate guess from the audience that primes might be special. But that theory was quickly exploded by the case of multiplier 21.

Table[Labeled[   PadLeft[IntegerDigits[      NestList[a FromDigits[RotateLeft[IntegerDigits[#, 2]], 2] + 1 &,        1, 100], 2]] // ArrayPlot, a], {a, 20, 24}]

OK, so then the hunt was on for what was special about the multipliers that led to complex behavior. But first we had to figure out how to recognize complex behavior. I thought I’d try something newfangled: using machine learning to make a feature space plot of the images for different multipliers.

It was somewhat interesting—and a nice application of machine learning—but not immediately too useful. (To make it better, one would have to think harder about the feature extractor to use.)

imags = Table[    PadLeft[IntegerDigits[       NestList[a FromDigits[RotateLeft[IntegerDigits[#, 2]], 2] + 1 &,         1, 100], 2]] // Image, {a, 50}]; FeatureSpacePlot[imags]

So how could one tell from that which were the complex patterns? A histogram of entropies wasn’t obviously illuminating:

Histogram[Entropy /@ imags]

As I was writing this blog post, I thought I should find the entropy distribution more accurately; even including 1000 possible multipliers, it still doesn’t seem terribly helpful:

Histogram[  Entropy /@    Table[PadLeft[      IntegerDigits[       NestList[a FromDigits[RotateLeft[IntegerDigits[#, 2]], 2] + 1 &,         1, 100], 2]] // Image, {a, 1000}]]

An expert in telecom math in the front row suggested taking a Fourier transform. I said I wasn’t hopeful:

Image[Abs[Fourier[ImageData[]]]]

Yes, there are better ways to do the Fourier transform. But someone else (a hedge-fund CEO, as it happened) suggested looking at the occurrences of particular 2×2 blocks in each pattern. For the case of multiplier 13, lots of blocks occur:

Counts[Flatten[   Partition[    PadLeft[IntegerDigits[      NestList[13 FromDigits[RotateLeft[IntegerDigits[#, 2]], 2] + 1 &,        1, 20], 2]], {2, 2}], 1]]

But for the case of multiplier 5, where the pattern is simple, most blocks never occur:

Counts[Flatten[   Partition[    PadLeft[IntegerDigits[      NestList[5 FromDigits[RotateLeft[IntegerDigits[#, 2]], 2] + 1 &,        1, 20], 2]], {2, 2}], 1]]

So this suggested that we just generate a list of how many of the 16 possible blocks actually do occur, for each multiplier:

blks = Table[   Length[Union[     Flatten[Partition[       PadLeft[IntegerDigits[         NestList[          a FromDigits[RotateLeft[IntegerDigits[#, 2]], 2] + 1 &, 1,           20], 2]], {2, 2}], 1]]], {a, 50}]

Here’s a plot:

ListLinePlot[blks]

Where are the 16s?

Flatten[Position[blks, 16]]

FindSequenceFunction didn’t have any luck with these numbers. Plotting the “block count” for longer gave this:

ListLinePlot[  Table[Length[    Union[Flatten[      Partition[       PadLeft[IntegerDigits[         NestList[          a FromDigits[RotateLeft[IntegerDigits[#, 2]], 2] + 1 &, 1,           20], 2]], {2, 2}], 1]]], {a, 1000}]]

Definitely some structure. But it’s not clear what it is.

And once again, we were out of time—having found an interesting kind of system with the curious property that it’s usually complex in its behavior, but for some special cases, isn’t.

The Live Experiment Process

I’ve done many live experiments over the years—though it’s been a while since they were about math. And as the Wolfram Language has evolved, it’s become easier and easier to do the experiments nicely and smoothly—without time wasted on glitches and debugging.

Wolfram Notebooks have the nice little feature that they (by default) keep a Notebook History (see the Cell menu)—that shows when each cell in the notebook has been modified. Here are the results for Experiment #1 and Experiment #2. Mostly they show rather linear progress, with comparatively little backtracking. (There’s a gap in Experiment #2, which came because my network connection suddenly stopped working. Conveniently, there were some networking experts in the audience—and eventually it was determined that the USB-C connection from my fine new computer to the projector had somehow misnegotiated itself as an Ethernet connection…)

Cell history

Every year at our Summer School I start out by doing a live experiment or two—because I think live experiments are a great way to show just how accessible discovery can be if one approaches it the right way, and with the right tools. I’m expecting that live experiments will be an important part of the process of educating people about computational thinking too.

With the Wolfram Language, one can do live experiments—and live coding—about all sorts of things. (We even tried doing a prototype Live Coding Competition at our last Wolfram Technology Conference; it worked well, and we’ll probably develop it into a whole program.)

But whether they’re live or not, computer experiments are an incredibly powerful methodology for making discoveries—not least in mathematics.

Of course, it’s easy to generate all kinds of random facts about mathematics. The issue is: how does one generate “interesting” facts? In a first approximation, for a fact to be interesting to us humans, it has to relate to things we care about. Those things could be technological applications, observations about the real world—or just pieces of mathematics that have, for whatever reason, historically been studied (think Fermat’s Last Theorem, for example).

I like to think that my book A New Kind of Science significantly broadened the kinds of “math-like facts” that one might consider “interesting”—by providing a general intellectual framework (about computation, complexity, and so on) into which those facts can be fit.

But part of the skill needed to do good experimental mathematics is to look for facts that somehow can ultimately be related to larger frameworks, and ultimately to the traditions of mathematics. Like in any area of research, it takes experience and intuition—and luck can help too.

But in experimental mathematics, it’s extremely easy to get started: there’s plenty of fertile territory to be explored, even with quite elementary mathematical ideas. We just happen to live at a time when the tools to make this kind of exploration feasible first exist. (Of course, I’ve spent a lot of my life building them…)

How should experimental mathematics be done? Perhaps there could be “math-a-thons” (or “discover-a-thons”), analogous to hackathons, where the output is math papers, not software projects.

More than 30 years ago I started the journal Complex Systems—and one of my long-term goals was to make it a repository for results in experimental mathematics. It certainly has published plenty of them, but the standard form of modern academic papers isn’t optimized for experimental mathematics. Instead, one can imagine some kind of “Discoveries in Experimental Mathematics,” that is much more oriented towards straightforward reports of the results of experiments.

In some ways it would be a return to an earlier style of scientific publishing—like all those papers from the 1800s reporting sighting of strange new animals or physical phenomena. But what’s new today is that with the Wolfram Language—and particularly with Notebook documents—it’s possible not just to report on what one’s seen, but instead to give a completely reproducible version of it, that anyone else can also explore. (And if there’s a giant computation involved, to store the results in a cloud object.)

I’m hoping that finally it’ll now be possible to establish a serious ecosystem for experimental mathematics. A place where results can be found, presented clearly with good visualization and so on, and published in a form where others can build on them. It’s been a long time coming, but I think it’s going to be an important direction for mathematics going forward.

And it was fun for me (and I hope for the audience too) to spend a couple of hours prototyping it live and in public a few days ago.


Download the complete notebooks:
Session #1/Experiment #1 »
Session #2/Experiment #2 »

Try the experiments for yourself in the Wolfram Open Cloud »

1 comment

  1. Thank you so much for sharing the ideas behind your MoMath talk. I loved everything about this post and used it for a fun project with my kids tonight:

    https://mikesmathpage.wordpress.com/2017/03/07/sharing-stephen-wolframs-momath-talk-with-kids/