Thursday, November 15, 2012

Welcome!

Hey everyone! This is my second blog post. I just finished up my first ever blog post about an hour ago. I've never blogged before. In the back of my mind, it's always seemed a little self-indulgent to me to think that everyone should care enough about whatever I'm thinking to come read about it. Nevertheless, here I am, a "blogger". 

Who knows what all I will end up using this blog for, but right now I'm planning on using it to explain my random musings about math, programming, and science. That might sound dry and boring to you. So my goal is to have at least one simple statement or image in every post that absolutely anyone can look at and go "hey, that's pretty cool." If my posts do that, then this blog will be a success to me.

I'm going to try to keep the posts light on jargon. So I might not always use the "technically correct" terminology or notation, merely the one that I think is most agreeable and understandable.

If you want to get a flavor of the kinds of things I'm going to be posting about, scroll down and check out my first post! Hopefully you'll go, "hey, that's pretty cool." Maybe you'll even learn something.

The Collatz Conjecture: emergent order in numbers

Hey everyone, welcome to my first blog post. This blog post is all about the beautiful emergent order in sequences of numbers from the Collatz conjecture. Not familiar with the Collatz conjecture? Neither was I until today. Read on. If you're not at all interested in a tiny math lesson, then scroll down to the images. There's a couple cool gifs that look kind of pretty (in my opinion). But what really makes them cool is how we arrive at them, so read on.

I was hitting "random" on my favorite comic strip website xkcd.com when I stumbled across this one: xkcd.com/710/. I was unfamiliar with the "Collatz conjecture", so I went to the internet's #1 resource on quick information referencing to learn more about it (in case this is your first time on the internet, I'm referring to Wikipedia: en.wikipedia.org/wiki/Collatz_conjecture).

The Collatz conjecture is one of those quirky numerical observations that seems completely random at first glance, but has managed to spark some mathematical interest because it allows us to probe our understanding of the very foundations of mathematics - sequences of numbers themselves. Paraphrased, the Collatz conjecture states the following: 

Start with any positive integer n. If n is even, divide it in half. If it is odd, multiply it by three and add one. If you continue this process for a sufficient number of steps, you will eventually obtain n=1, regardless of what n started as.

The conjecture is simple enough; a third grader should be able to carry out this procedure on a given number. However, there is no obvious reason why this conjecture should necessarily be true or false, and indeed (according to Wikipedia), the problem of whether or not the conjecture is true is "algorithmically undecidable".

Let's take a look at a few simple examples of the Collatz conjecture in action. 

If we start with 4, we see that n quickly converges to 1. The starting number 4 is divided in half to become 2, which is divided in half to become 1. I will refer to this as the "trajectory" of 4, which we can represent as {4,2,1}. If we start with 3, the trajectory is {3,10,5,16,8,4,2,1}. 

I will refer to the number of steps it takes n to converge to 1 as the "trajectory length". So, 4 has a trajectory length of 2 and 3 has a trajectory length of 7. 

If we plot starting numbers n vs trajectory lengths for all numbers between, say, 1 and 10, we don't see any remarkable pattern: 



If we plot trajectory lengths for starting numbers between 1 and 100, some sort of pattern starts to emerge, though it is perhaps not what we would expect.


Let's stop with the cute stuff: let's plot all the way up to 1 million:


Woah. If this graph doesn't hit you over the head with the beauty that has emerged from its randomness, then you are reading the wrong blog.

To illustrate the emergent order better, I made a gif image of the trajectory lengths graphs, with k increasing.



Are there any other emergent patterns in the Collatz conjecture? Definitely. Remember that every starting number n has a trajectory, and for instance the trajectory of 3 is {3,10,5,16,8,4,2,1}. The max of the trajectory I will define as the maximum value that the trajectory crosses. So the max of the trajectory of 3 is 16. Let's plot the max of the first 50 Collatz trajectories.

So far, we don't see anything too remarkable. The max tends to increase with larger starting values, as we might expect, but there's no remarkable trend. Let's cut to the chase and plot the max of the Collatz trajectories of the first 1 million numbers.


Again, we see an amazing an unexpected order emerge. Numbers are amazing.

For your viewing pleasure, I again made a gif image out of this data so that you can watch the pattern emerge and grow:





If you're not interested in programming, you can probably stop reading here. The rest of this post is going to indulge my inner programmer (and Python lover).

So how did I make these graphs? With the aid of my favorite programming tool, PYTHON! 

Wikipedia gives us some simple pseudocode to trace the trajectories (which is really just a restatement of the conjecture itself). This can be easily translated into a Python function which returns the length of the trajectory of n. (As a side note, this function assumes that the conjecture is true for all input numbers! If it's not true for some input number, then we'll get an infinite while loop. A more robust function would test if the trajectory has formed a closed cycle, though this would probably slow the function down quite a bit).

def traj(n):
    ops = 0
    while n > 1:
        ops = ops + 1
        if n % 2 == 0:
            n=n/2
        else:
            n=3*n+1
    return ops

We can embed this function in a for loop that stores the trajectory lengths between for starting numbers between 1 and k:

results = [0,0]
for i in range(2,n):
        results.append(traj(i))

and plot them using the Python package Matplotlib (if you happen to be following along, you will need the Matplotlib package installed):

import matplotlib.pyplot as plots
ms=2
plots.plot(results,'o',ms=ms)
plots.ylabel("trajectory length")
plots.xlabel("starting number, n")
plots.show()


How far out can we plot using my simple Python code? Well, it depends how long you'd like to wait for it to run. We can time how long it takes our code to run using time.clock() in Python. If we time our code for various k (the max integer whose trajectory length we are computing), and plot on a log scale, we get the following graph:


Computing a million trajectories took about 20 seconds. We can see that for every factor of 10 that we increase the number of trajectories by, we increase the computing time by about a factor of 10. At k=1,000,000, it takes about 20 seconds to run. At k=1,000,000,000, we could estimate that it would take about 20,000 seconds to run (several hours). I don't know about you, but I don't really like waiting that long for my code to run. Can we speed up our code at all? You bet we can!

One algorithmic simplification we can make is instead of calculating 3*n+1, we can calculate (3n+1)/2, since we know that 3n+1 will always be greater than 1 anyway. We just need to make sure we record (3n+1)/2 as 2 operations. So our new traj function would look like:

def traj(n):
    ops = 0
    while n > 1:
        if n % 2 == 0:
            n=n/2
            ops=ops+1
        else:
            n=(3*n+1)/2
            ops=ops+2
    return ops

Does this speed up our computation? Not much, but a bit. In this graph, the blue points are the computation times without our new optimization, and the green points are the computation times with our new optimization. We've speeded things up a bit, but I still want to get faster.


I figured out a more substantial optimization when I was looking at the image that inspired this whole inquiry (http://xkcd.com/710/). The fact is, a lot of trajectories meet up, and once they meet each other, they never leave each other (thus a graph representation for the trajectories is convenient). For instance, the trajectory for n=6 is {6,3,10,5,16,8,4,2,1}, while the trajectory for n=3 is {3,10,5,16,8,4,2,1}. So if we've already computed the length of the trajectory of 3, then we can quickly figure out that the trajectory length of 6 is only one more than the trajectory length of 3 because the first operation in the trajectory of n=6 is n=6/2=3. The code that would use this information to our advantage would look like this (in basic machine learning, this is called "memorization", I believe):

results = [0,0]
for i in range(2,n):
    results.append(trajLearner(i,results))



def trajLearner(n, learner):
    init = n
    ops = 0
    while n > 1:
        if n < init:
            ops = ops + learner[n]
            break
        else:
            if n % 2 == 0:
                n=n/2
                ops=ops+1
            else:
                n=(3*n+1)/2
                ops=ops+2
    return ops

This optimization improves speed, on average, by about a factor of 10. The red dots are computation times with our new optimization. 



There are more optimizations we could do (and I may do a little later. If we really cared about speed, I'd code it up in C, but I like Python too much), but I think this is more than sufficient for my first blog post. I hope you enjoyed reading my first post as much as I enjoyed writing it! I'll probably add more to this post when I get around to it. Let me know if you have any feedback.

Mike