A Neighborhood of Infinity

Some parallels between classical and quantum mechanics

2012-01-21T13:38:00.000-08:00

Introduction
This isn't really a blog post. More of something I wanted to interject in a discussion on Google plus but wouldn't fit in the text box.

I've always had trouble with the way the Legendre transform is introduced in classical mechanics. I know I'm not the only one. Many mathematicians and physicists have recognised that it seems to be plucked out of a hat like a rabbit and have even written papers to address this issue. But however much an author attempts to make it seem natural, it still looks like a rabbit to me.

So I have to ask myself, what would make me feel comfortable with the Legendre transform?

The Legendre transform is an analogue of the Fourier transform that uses a different semiring to the usual. I wrote briefly about this many years ago. So if we could write classical mechanics in a form that is analogous to another problem where I'd use a Fourier transform, I'd be happier. This is my attempt to do that.

When I wrote about Fourier transforms a little while back the intention was to immediately follow it with an analogous article about Legendre transforms. Unfortunately that's been postponed so I'm going to just assume you know that Legendre transforms can be used to compute inf-convolutions. I'll state clearly what that means below, but I won't show any detail on the analogy with Fourier transforms.

Free classical particles
Let's work in one dimension with a particle of mass whose position at time is . The kinetic energy of this particle is given by . Its Lagrangian is therefore .

The action of our particle for the time from to is therefore

The particle motion is that which minimises the action.

Suppose the position of the particle at time is and the position at time is . Then write for the action minimising path from to . So

where we're minimising over all paths such that .

Now suppose our system evolves from time to . We can consider this to be two stages, one from to followed by one from to . Let be the minimised action analogous to for the period to . The action from to is the sum of the actions for the two subperiods. So the minimum total action for the period to is given by

Let me simply that a little. I'll use where I previously used and for . So that last equation becomes:

Now suppose is translation-independent in the sense that . So we can write . Then the minimum total action is given by

Infimal convolution is defined by

so the minimum we seek is

So now it's natural to use the Legendre transform. We have the inf-convolution theorem:

where is the Legendre transform of given by

and so (where we use to represent Legendre transform with respect to the spatial variable).

Let's consider the case where from onwards the particle motion is free, so . In this case we clearly have translation-invariance and so the time evolution is given by repeated inf-convolution with and in the "Legendre domain" this is nothing other than repeated addition of .

Let's take a look at . We know that if a particle travels freely from to over the period from to then it must have followed the minimum action path and we know, from basic mechanics, this is the path with constant velocity. So

and hence the action is given by

So the time evolution of is given by repeated inf-convolution with a quadratic function. The time evolution of is therefore given by repeated addition of the Legendre transform of a quadratic function. It's not hard to prove that the Legendre transform of a quadratic function is also quadratic. In fact:

Addition is easier to work with than inf-convolution so if we wish to understand the time evolution of the action function it's natural to work with this Legendre transformed function.

So that's it for classical mechanics in this post. I've tried to look at the evolution of a classical system in a way that makes the Legendre transform natural.

Free quantum particles
Now I want to take a look at the evolution of a free quantum particle to show how similar it is to what I wrote above. In this case we have the Schrödinger equation

Let's suppose that from time onwards the particle is free so . Then we have

Now let's take the Fourier transform in the spatial variable. We get:

We can write this as

where

So the time evolution of the free quantum particle is given by repeated convolution with a Gaussian function which in the Fourier domain is repeated multiplication by a Gaussian. The classical section above is nothing but a tropical version of this section.

Conclusion
I doubt I've said anything original here. Classical mechanics is well known to be the limit of quantum mechanics as and it's well known that in this limit we find that occurrences of the semiring are replaced by the semiring . But I've never seen an article that attempts to describe classical mechanics in terms of repeated inf-convolution even though this is close to Hamilton's formulation and I've never seen an article that shows the parallel with the Schrödinger equation in this way. I'm hoping someone will now be able to say to me "I've seen that before" and post a relevant link below.

Note
I'm not sure how the above applies for a non-trivial potential . I wrote this little Schrödinger equation solver a while back. As might be expected, it's inconvenient to use the Fourier domain to deal with the part of the evolution due to . In order to simulate a time step of the code simulates in the Fourier domain assuming the particle is free and then spends solving for the -dependent part in the spatial domain. So even in the presence of non-trivial it can still be useful to work with a Fourier transform. Almost the same iteration could be used to numerically compute the action for the classical case.

Lossless decompression and the generation of random samples

2012-01-07T09:40:00.000-08:00

Let S be some finite set with a probability distribution on it. Here's a diagram showing some example probabilities for the set S={A, B, C, D, E}:

How can we generate lots of random samples from this distribution? A popular method involves first storing the cumulative distribution function in a table like so:

We then use any of a number of popular methods to generate uniform pseudorandom numbers in the range [0,1) and for each one walk through the table until we find the first entry greater than our pseudorandom number. The symbol above the number is the one we generate. So if we generated 0.512, we'd pick symbol C. It's straightforward to prove this gives the correct probability distribution.

As described this algorithm can be slow. If the size of the table is N we may have to walk through up to N entries to generate each sample.

One approach to accelerating this algorithm is to quantise our pseudorandom number and use it to look up, in a precomputed table, a jump into our cumulative distribution table. I've used this several times in my visual effects career. But today I'm going to take a different approach.

Another natural way to speed up sample generation is to use a binary search to find the appropriate point in the cumulative distribution, for example the C++ standard template library's upper_bound method will do the job.

But what is a binary search algorithm going to do? Typically it's going to start by comparing our pseudorandom number with the midpoint of the table. If our number is bigger then it will recursively use binary search on the left (looking next at the midpoint of the left half), otherwise on the right, and so on. If we're generating many samples from the same distribution we're going to be repeatedly looking up the same midpoints in the table. At the end of the day, the process can be described by a decision tree. So we may as well throw away the table and build the decision tree up front. Here's what it might look look like:

But maybe we can do better. For example C is three times as likely as A but they both take the same amount of time to generate as they both require walking down to depth 2. Meanwhile D has a probability of 0.25 and requires walking to depth 3. We'd like to rebalance the tree. Note also that there's no reason to list our original PDF in the order I gave. Different orderings might give different trees.

It's straightforward to describe what we want to optimise. We want to place sample i at depth d_i so that the expected value of d_i is as small as possible. In other words we want to minimise Σp_id_i. But this is precisely the problem solved by Huffman coding.

And that's the conclusion I wanted to reach: Huffman coding gives the optimal decision tree to generate random samples. It also tells us this interesting consequence: if we use a decision tree method then the performance of our algorithm is bounded by the entropy of the probability distribution. I find this connection between entropy and algorithmic complexity pretty surprising.

I learnt the above during my interview at Google!

Why is there this connection between a compression algorithm and the generation of random samples? It took me a little longer to realise why but it's quite straightforward.

Huffman coding tries to compress text one letter at a time on the assumption that each letter comes from some fixed and known probability distribution. If the algorithm is successful then we'd expect the compressed text to look like a uniformly distributed sequence of bits. If it didn't then there'd be patterns that could be used for further compression.

So when we're decompressing Huffman encoded data we have a machine that takes as input uniformly distributed bits and which outputs letters sampled from some probability distribution. But that's exactly the same problem that I posed above. So, at least from some perspective, decompression is precisely the same thing as generating samples from a probability distribution.

Or to put it another way: there is a class of algorithm whose purpose is to convert samples from one probability distribution into samples from another. When we use one of these algorithms to convert from samples from some distribution P that's easy to generate samples from, into samples from Q, we call this "an algorithm to sample from distribution Q". When we use one of these algorithms to convert from some distribution to a uniform distribution, and the function given by the algorithm is invertible, then we call it "lossless compression".

Quick and dirty reinversion of control

2011-10-30T13:19:00.000-07:00

It's taken for granted by many people that Haskell and static types are incompatible with prototyping and quick-and-dirty hacks.

I wanted to put together some OpenGL code that had a script for how a bunch of graphics should be displayed. It was essentially an imperative specfification for my program. For quick and dirty hacks, GLUT is tolerable. But even when programming in C/C++, it's not supportive of programming in a straightforward imperative style because it uses inversion of control. Many graphics applications are written in a state machine style where the state machine gets to tick once each time there is an event callback. This really doesn't fit the imperative script style.

But it is possible to reinvert inversion of control in any language that supports continuations. And that includes languages like Python that support linear continuations in the form of generators. But I'm using Haskell here.

Continuations reify the remainder of a computation. Or in more down to earth language: they allow you to grab the stuff you're about to do as a function, put it on ice for a while, and then carry on doing it later. So imagine we had a block of imperative code and that we'd like, at each GLUT callback, to make some progress through this block. We can use continuations like this: each time we want to yield control back to the main loop we simply grab the remainder or our 'script' as a continuation and make it the callback to be executed next time GLUT is ready.

The slight wrinkle is that OpenGL/GLUT calls use IO. To combine IO and continuations we need the ContT monad transformer.

I'll do everything except the yield function first and get back to that at the end.

Some standard library imports:

> import Graphics.UI.GLUT
> import Control.Monad.Cont

Some simple code to draw a line from left to right:

> display :: GLdouble -> IO ()
> display y = do
>   clear [ColorBuffer]
>  
>   renderPrimitive LineStrip $ do
>     vertex (Vertex2 (-1) (-y))
>     vertex (Vertex2 1 y)

>   swapBuffers
>   postRedisplay Nothing

Some standard OpenGL/GLUT setup:

> main = do
>   (progname, _) <- getArgsAndInitialize
>   initialWindowSize  $= Size 500 500
>   initialDisplayMode $= [DoubleBuffered, RGBMode]
>   createWindow "Bounce!"

>   matrixMode $= Modelview 0
>   loadIdentity

>   matrixMode $= Projection
>   loadIdentity
>   ortho (-1) 1 (-1) 1 (-1) 1

Our script is called before the main loop.

>   imperative
>   mainLoop

And now comes the actual script. Apart from the liftIO calls this should be almost as easy to read as BASIC programming from the days of yore:

> imperative = flip runContT return $ do

>   liftIO $ print "Start!"

>   forever $ do

>     forM_ [-1, -0.992 .. 1.0] $ \y -> do
>       render $ display y
>       yield

>     liftIO $ print "Bounce!"

>     forM_ [-1, -0.992 .. 1.0] $ \y -> do
>       render $ display (-y)
>       yield

>     liftIO $ print "Bounce!"
>     yield

The first thing to note is that render doesn't actually do any rendering. At the end of the day we can't tell GLUT when to render, it only calls you. So instead render tells GLUT what to do next time it's in the mood for a bit of rendering:

> render f = liftIO $ displayCallback $= f

That leaves one thing to explain: yield. It needs to grab the remainder of the script and package it up in a form suitable for installation as an idle callback. But there's a catch: continuations are notorious for making your head explode. If you're throwing together a quick and dirty hack, that's the last thing you need. Here's where static types come to the rescue. As Conor McBride points out, we want to just do the easy thing and follow gravity downwards.

So first we try to guess the type of yield. We know we're working with the ContT IO monad. So its type is going to be ContT IO a for some a. There's no particular type of data we want to get out of yield, it's just a thing we want executed. So we can guess the type is ContT IO (), the type () being the generic filler type when we don't actually have any data.

Let's look at the definition of ContT:

newtype ContT r m a = Cont {
    runContT :: (a -> m r) -> m r
}

The type r is the final return type from our continuation. We're not interested in a return value, we just want to *do* stuff. So we expect r to be () as well.

So yield must essentially be of type (() -> IO ()) -> IO ().

So we want to concoct something of this type using GLUT's idleCallback function. As yield must take a function as argument it must look something like:

yield = ContT $ \f -> ...

We know that f is of type () -> IO (). So there's only one thing we can do with it: apply it to (). That gives us something of type IO (). That's precisely the type of GLUT's idleCallback. So we put it all together:

> yield = ContT $ \f -> idleCallback $= Just (f () )

The code now works. I didn't have to spend even a moment thinking about the meaning of a continuation. Implementing yield was about as hard as putting together a jigsaw puzzle with three pieces. There's only so many ways you can put the pieces together.

And that's a simple example of why I often like to write quick-and-dirty code with a statically typed language.

(Oh, and I'm not trying to take part in a war. I like to prototype in many different languages, some of which are dynamically typed.)

PS Note also that the above code illustrates one way to avoid IORefs in GLUT code.

Computing errors with square roots of infinitesimals.

2011-08-13T08:12:00.000-07:00

Abstract
Automatic differentiation (AD) gives a way to carry out uncertainty propagation. But used in the obvious way it leads to bias. This article introduces "square roots of infinitesimals" that can be used to give more accurate results.

Introduction
In the real world measurements have errors and we often want to know how much our final answers are affected by those errors. One tool for measuring the sensitivity to errors of our results is calculus. In fact, we can use automatic differentiation to give us a nice way to model error. Here's an implementation:

> import Control.Monad

> import Control.Monad.State

> import Control.Applicative

> import qualified Data.IntMap as I



> infixl 6 .+.

> infixl 7 .*

> infixl 7 *.



> data D a = D a a deriving (Eq, Show)



> instance Num a => Num (D a) where

>   D x a+D x' a' = D (x + x') (a + a')

>   D x a*D x' a' = D (x*x') (a*x' + x*a')

>   negate (D x a) = D (negate x) (negate a)

>   fromInteger n = D (fromInteger n) 0

>   abs _ = error "No abs"

>   signum _ = error "No signum"



> d = D 0 1

As I've talked about before, the value d can be thought of as an infinitesimal number whose square is zero. However, to first order we can replace d with a small number and use it to compute errors. Here's a function to perform such a substitution:

> approx :: Num a => D a -> a -> a

> approx (D x d) e = x+d*e

Suppose we have a square whose side we've measured as 1m to an accuracy of 1cm. We can represent this as:

> sq_side = D 1 0.01

We can now compute the area:

> sq_area = sq_side^2

We get D 1.0 2.0e-2. We can interpret this as meaning the area is 1m² with an accuracy of 0.02m².

We can make "to an accuracy of" more precise. Differentiation models a function locally as an approximate linear function. (1m+δ)²=1m²+2mδ+O(δ²). So at 1m, the function to compute area locally scales lengths by 2m. So if the length measurement is actually a sample from a distribution with given mean 1m and standard deviation 0.01m, the area is approximately a sample from a distribution with mean 1m² and SD 0.02m².

This approach is nice, but sometimes we want a little more accuracy in our estimates. In particular, if you square a sample from a normal distribution with small variance and positive mean, then the nonlinearity of the squaring operation means that samples that are larger than the mean move further away from the mean, when squared, than samples less than the mean. So we should actually expect our area computations to be slightly biased upwards from 1m². Unfortunately, this is a second order effect that isn't visible from looking only at first derivatives.

That's not a problem, we can easily compute second derivatives using automatic differentiation. However, that can complicate things. What happens if we use multiple measurements to compute a quantity? Each one is a different sample from a different distribution and we don't want these measurements to be correlated. If we approach this in the obvious way, when we want to use n measurements we'll need to compute n² partial second derivatives. However, by tweaking AD slightly we'll only need n derivatives.

Square roots of infinitesimals
In addition to the usual infinitesimal d we want to introduce quantities, w_i, that represent independent random "noise" variables that are infinitesimal in size. We'll be interested in expectation values so we'll also need an expectation function, e. We want e(w_i)=0. But w_i² is always positive so its expectation is always greater than or equal to zero. We want our random variables to be infinitesimal so we pick e(w_i²)=d. We also want e(w_iw_j)=0 because of independence. If the w_i are already infinitesimal, the dw_i should be zero. So let's define an algebraic structure that captures these relationships. So we extend D a by introducing w i so that:

Any element of this algebra can be written as x+ad+Σb_iz_i. We represent b sparsely by using an IntMap. Here's an implementation:

> data S a = S a a (I.IntMap a) deriving (Eq, Show)



> (.+.) :: Num a => I.IntMap a -> I.IntMap a -> I.IntMap a

> ( *.) :: Num a => a -> I.IntMap a -> I.IntMap a

> (.*)  :: Num a => I.IntMap a -> a -> I.IntMap a

> (.*.) :: Num a => I.IntMap a -> I.IntMap a -> a



> (.+.) = I.unionWith (+)

> a *. v = I.map (a *) v

> v .* b = I.map (* b) v

> a .*. b = I.fold (+) 0 $ I.intersectionWith (*) a b



> instance Num a => Num (S a) where

>   S x a b+S x' a' b' = S (x + x') (a + a') (b .+. b')

>   S x a b*S x' a' b' = S (x*x') (a*x' + x*a' + b.*.b') (x*.b' .+. b.*x')

>   negate (S x a b) = S (negate x) (negate a) (I.map negate b)

>   fromInteger n = S (fromInteger n) 0 I.empty

>   abs _ = error "No abs"

>   signum _ = error "No signum"

Here are the individual w i:

> w :: Num a => Int -> S a

> w i = S 0 0 (I.fromList [(i, 1)])

We compute expectation values linearly by mapping the w i to zero:

> e :: Num a => S a -> D a

> e (S x a _) = D x a

We can also represent numbers whose values we know precisely:

> sure x = S x 0 I.empty

Example
Let's revisit the area example. This time we can represent the length of the side of our square as

> sq_side' = 1+0.01*w 0

> sq_area' = sq_side'^2

We get S 1.0 1.0e-4 (fromList [(0,2.0e-2)]). We can directly read off that we have a bias of 10^-4m² which is 1cm^2. We can encapsulate this as:

> mean f = approx (e f) 1

We can directly read off the variance from the element of the algebra. However, we can also compute the variance using mean. It's just:

> var f = mean ((f-sure (mean f))^2)

(Note that this gives a very slightly different result from the value you can read off directly from the S object. It depends on whether we're measuring the deviation around the unbiased or biased mean. To the order we're considering here the difference is small. Here's var' anyway:

> var' (S _ _ v) = I.fold (+) 0 $ I.map (\x -> x^2) v

)

We can also define covariance:

> cov f g = mean ((f-sure (mean f))*(g-sure (mean g)))

More functions
We can now follow through just like with automatic differentiation to compute lots more functions. We use the fact that:

> instance Fractional a => Fractional (S a) where

>   fromRational x = S (fromRational x) 0 I.empty

>   recip (S x a b) = let r = recip x

>                     in S r (-a*r*r+r*r*r*(b.*.b)) ((-r*r)*.b)



> instance Floating a => Floating (S a) where

>   pi = sure pi

>   sin (S x a b) = let s = sin x

>                       c = cos x

>                   in S s (a*c - s/2*(b.*.b)) (c*.b)

>   cos (S x a b) = let s = sin x

>                       c = cos x

>                   in S c (-a*s - c/2*(b.*.b)) ((-s)*.b)

>   exp (S x a b) = let e = exp x

>                   in S e (a*e + e/2*(b.*.b)) (e*.b)

>   sqrt (S x a b) = let s = sqrt x

>                   in S s (a/(2*s)-1/(4*s*s*s)*(b.*.b)) (1/(2*s)*.b)

>   log (S x a b) = let r = 1/x

>                   in S (log x) (r*a-r*r/2*(b.*.b)) (r*.b)

>   asin = undefined

>   acos = undefined

>   atan = undefined

>   sinh = undefined

>   cosh = undefined

>   tanh = undefined

>   asinh = undefined

>   acosh = undefined

>   atanh = undefined

A real example
Let's make this effort worthwhile. We'll compute errors for a computation that uses the errors in a messy nonlinear way. Suppose we're in the lab measuring radioactive decay. We measure the geiger counter reading at times t = 0hr, 1hr, 2hr, 3hr, 4hr at which point we compute an estimate for when the decay will drop to one tenth of its original value. We'll assume the decay fits a model counts/sec = a exp(-λt) and that the counts have an error with SD 0.05. We're going to compute the error in the estimated time to hit one tenth radioactivity in the case when the half life is 30 minutes and a=2:

> t = [0..4]

> counts = map (\i-> 2*exp(-0.5*fromIntegral i)+0.05*w i) t

We'll be fitting a curve using logarithmic regression so we'll need the following function. Given a pair of lists x and y it returns (m, c) where y=mx+c is the standard least squares fit.

> regress :: Fractional a => [a] -> [a] -> (a, a)

> regress x y =

>     let sx = sum x

>         sy = sum y

>         sxx = sum $ map (^2) x

>         sxy = sum $ zipWith (*) x y

>         n = fromIntegral (length x)

>         s = 1/(sx*sx-n*sxx)

>     in (s*sx*sy-s*n*sxy, -s*sxx*sy+s*sx*sxy)

Logarithmic regression:

> (m, c) = regress (map fromIntegral t) (map log counts)

> lambda = -m

> a = exp c

> t_tenth = -log (0.1/a)/lambda

We can now go ahead and compute the mean and variance of our estimate:

*Main> mean t_tenth 

5.98036172868899

*Main> var t_tenth

0.15583537298560224

The correct time is about 5.991 so the regression method above is biased by about 0.01. If we repeated the same experiment over and over again and averaged the estimates we got from logarithmic regression the process would not converge to the correct result. In fact, we can compute "ground truth" by simulating the experiment a million times in Octave and estimate the mean and variance from that. The code is in the appendix. Obviously this is a much slower process but it clearly demonstrates the biasedness of using regression this way.

GNU Octave, version 3.4.0

Copyright (C) 2011 John W. Eaton and others.

ans =  5.9798

ans =  0.15948

Final thoughts
This is yet another example of extending automatic differentiation. We have variants for single variable differentiation, multivariate differentiation, multiple differentiation, divided differences, splitting a function into odd and even parts and now automatic error propagation.

This stuff was very loosely inspired by reading An Introduction to Stochastic Processes in Physics. I'm attempting to capture the semi-formal rules used in that book to reason about differentials and you can think of the algebra above as representing stochastic differentials. I made a guess that the algebra is called the Itō algebra. Sure enough, you'll get a few hits.

The most similar published work I can find is Automatic Propagation of Uncertainties but it seems to just use ordinary AD.

This technique may be useful for Extended Kalman Filtering.

I haven't done the work to make precise statements about how accurate you can expect my estimates of expectations to be.

It's possible to implement a monad with syntax similar to other probability monads by using state to bump up the i in w i each time you generate a new random variable. But bear in mind, these are always intended to be used as *infinitesimal* random variables.

Appendix: Octave code

m = 5;

n = 1000000;



x = repmat([0:m-1]',1,n);

y = repmat([2*exp(-0.5*[0:m-1]')],1,n)+0.05*normrnd(0,1,m,n);



sx = sum(x);

sxx = sum(x.*x);

p = sum(log(y));

q = sum(x.*log(y));



s = 1./(sx.*sx-m*sxx);

m = s.*sx.*p-m*s.*q; # Redefined

c = -s.*sxx.*p+s.*sx.*q;



lambda = -m;

a = exp(c);

x_tenth = -log(0.1./a)./lambda;



mean(x_tenth)

var(x_tenth)

Profunctors in Haskell

2011-07-23T14:39:00.000-07:00

> {-# LANGUAGE TypeSynonymInstances, RankNTypes, ExistentialQuantification #-}

Introduction
When I wrote about coends a while back I made up a term 'difunctor'. More recently it was pointed out to me that the correct word for this concept is 'profunctor', but unfortunately my knowledge came from MacLane which mentions that word nowhere.

Profunctors are ubiquitous in Haskell programming. Probably the most natural definition of Hughes Arrows is via profunctors. Profunctors also play a role a little like tensors leading to a use of the terms 'covariant' and 'contravariant' that looks remarkably like the way those terms are used in tensor calculus.

For categories C and D, A profunctor is a functor D^op×C→Set and is written C↛D. (I hope that arrow between C and D is in your font. It's missing on iOS.)

I'll reuse my Haskell approximation to that definition:

> class Profunctor h where
>   lmap :: (d' -> d) -> h d c -> h d' c
>   rmap :: (c -> c') -> h d c -> h d c'

We need cofunctoriality for the first argument and functoriality for the second:

lmap (f . g) == lmap g . lmap f
rmap (f . g) == rmap f . rmap g

(Strictly we probably ought to call these 'endoprofunctors' as we're only really dealing with the category of Haskell types and functions.)

There are lots of analogies for thinking about profunctors. For example, some people think of them as generalising functors in the same way that relations generalise functions. More specifically, given a function f:A→B, f associates to each element of A, a single element of B. But if we want f to associate elements of A with elements of B more freely, for example 'mapping' elements of A to multiple elements of B then we instead use a relation which can be written as a function f:A×B→{0,1} where we say xfy iff f(x,y)=1. In this case, profunctors map to Set rather than {0,1}.

A good example is the type constructor (->)

> instance Profunctor (->) where
>   lmap f g = g . f
>   rmap f g = f . g

It's common that the first argument of a profunctor describes how an element related to a type is sucked in, and the second describes what is spit out. a -> b sucks in an a and spits out a b.

Given a function f we can turn it into a relation by saying that xfy iff y=f(x). Similarly we can turn a functor into a profunctor. Given a functor F:C→D we can define a profunctor F^*:C↛D by

> data UpStar f d c = UpStar (d -> f c)
> instance Functor f => Profunctor (UpStar f) where
>   lmap k (UpStar f) = UpStar (f . k)
>   rmap k (UpStar f) = UpStar (fmap k . f)

You may be able to see how the second argument to a profunctor sort of plays a similar role to the return value of a functor, just as the second argument to a relation sometimes plays a rule similar to the return value of a function.

There also an opoosing way to make a profunctor from a functor just as there is with functions and relations:

> data DownStar f d c = DownStar (f d -> c)
> instance Functor f => Profunctor (DownStar f) where
>   lmap k (DownStar f) = DownStar (f . fmap k)
>   rmap k (DownStar f) = DownStar (k . f)

Note that the identity functor gives us something isomorphic to (->) whether you use UpStar or DownStar.

Dinatural transformations
Just as we have natural transformations between functors, we have dinatural transformations between profunctors. My previous definition of dinatural was specialised to a particular case - dinaturals between a profunctor and the constant profunctor.

Firstly, let's think about natural transformations. If F and G are functors, and h is a natural transformation h:F⇒G, then we have that

h . fmap f = fmap f . h

If we think of F and G as containers, then this rule says that a natural transformation relates the structures of the containers, not the contents. So using f to replace the elements with other elements should be invisible to h and hence commute with it.

Something similar happens with dinatural transformations. But this time, instead of relating the argument to a natural transformation to its return result, it instead relates the two arguments to a profunctor.

Given two profunctors, F and G, A dinatural transformation is a polymorphic function of type:

> type Dinatural f g = forall a. f a a -> g a a

but we also want something analogous to the case of natural transformations. We want to express the fact that if phi :: Dinatural F G, then phi doesn't see the elements of F a a or G a a. Here's a way to achieve this. Suppose we have a dinatural transformation:

phi :: Dinatural G F

and a function f :: X -> X' then we can use lmap to apply f on the left or right of F and G. The definition of dinaturals demands that:

rmap f . phi . lmap f = lmap f . phi . rmap f

ie. that we can apply f on the left before applying phi, and then do f on the right, or vice versa, and still get the same result.

I'm not sure but I think that we don't need to check this condition and that just like the case of naturals it just comes as a free theorem.

Composing profunctors
It's easy to see how to compose functors. A functor is a polymorphic function from one type to another. It's not straightforward to compose profunctors. It's tempting to say that a profunctor maps a pair of types to a type so they can be composed like functions. But the original definition says the definition is D^op×C→Set. So as a function it doesn't map back to the category but to Set. For Haskell we replace Set with Hask, the category of Haskell functions and types. So we have Hask^op×Hask→Hask. It's easy invent a scheme to compose these because Hask appears 3 times. But it'd be wrong to exploit this in a general definition applying to many categories because in the proper definition of profunctor we can't assume that a profunctor maps back to the spaces you started with.

We can try composing profunctors by analogy with composing relations. Suppose R and S are relations. If T=S○R is the composition of R and S then xTz if and only if there exists a y such that xRy and ySz. If our relations are on finite sets then we can define T(x,z) = Σ_yR(x,y)S(y,z) where we work in the semiring on {0,1} with 0+0=0, 0+1=1+0=1+1=1 but with the usual product.

There is an analogue of "there exists" in Haskell - the existential type. Remembering that we write Haskell existential types using forall we can define:

> data Compose f g d c = forall a.Compose (f d a) (g a c)

As mentioned above, functors give rise to profunctors. It'd be good if composition of functors were compatible with composition of profunctors. So consider

Compose (UpStar F) (UpStar G)

for some F and G. This is essentially the same as

exists a. (d -> F a, a -> G c)

What can we discover about an element of such a type? It consists of a pair of functions (f, g), but we can't ever extract the individual functions because the type of a has been erased. To get anything meaningful out of g we need to apply it to an a, but we don't have one immediately to hand, after all, we can't even know what a is. But we do have an F a if we can make a d. So we can use fmap to apply g to the result of a. So we can construct fmap g . f :: d -> F (G c). There is no other information we can obtain. So the composition is isomorphic to UpStar of the functorial composition of F and G. Again, we can probably make this a rigorous proof by making use of free theorems, but I haven't figured that out yet.

But there's a catch: I said I wanted a definition that applies to more categories than just Hask. Well we can replace exists a with the coend operator. We also implicitly used the product operation in the constructor Compose so this definition will work in categories with suitable products. Symmetric monodial categories in fact.

Under composition of profunctors, (->) is the identity. At least up to isomorphism. This composition of profunctors is also associative up to isomorphism. Unfortunately the "up to isomorphism" means that we can't make a category out of profunctors in the obvious way. But we can make a bicategory - essentially a category where we have to explicitly track the isomorphisms between things that are equal in ordinary categories.

Profunctors as tensors
Given a profunctor F we can write F i j suggestively as F_i^j. Let's write the composition of F and G as ∃k. F_i^k G_k^j. We can use the Einstein summation convention to automatically 'contract' on pairs of upper and lower indices and write the composition as F_i^k G_k^j. The analogy is even more intriguing when we remember that in tensor notation, the upper indices are covariant indices and the lower ones are contravariant indices. In the case of profunctors, the two arguments act like the arguments to covariant and contravariant functors respectively. Note alse that because (->) is essentially the identity, we have →_i^jF_j^k=F_i^k. So (->) acts like the Kronecker delta. You can read more about this at mathoverflow where it is hinted that this analogy is not yet well understood. Note that we're naturally led to the trace of a profunctor: exists a. F a a.

Arrows as profunctors
The last thing I want to mention is that Hughes' Arrows are profunctors. There is an intuition that fits. If A is an Arrow, we often think of A d c as consuming something related to type d and emitting something related to type c. The same goes for profunctors. The full paper explaining this is Asada and Hasuo's Categorifying Computations into Components via Arrows as Profunctors with the profunctorial definition of Arrows given as Definition 3.2 (though that definition also appears in some earlier papers.)

The Infinitude of the Primes

2011-07-16T16:16:00.000-07:00

Introduction

> import Data.List hiding (intersect)

A colleague at work reminded me of Fürstenberg's topological proof of the infinitude of primes. But as I tried to argue a while back, topology is really a kind of logic. So Fürstenberg's proof should unpack into a sort of logical proof. In fact, I'm going to unpack it into what I'll call the PLT proof of the infinitude of the primes. I apologise in advance that I'm just going to present the unpacked proof, and not how I got there from Fürstenberg's.

A small formal language
We're going to start with a little language. Propositions of this language are of type Prop:

> data Prop = Modulo Integer Integer

The intention is that Modulo k n is the property of an integer being equal to k modulo n. More precisely, it represents the property of being writable in the form sn+k for some s. (We disallow n=0.) But I also want to allow people to combine properties using "and" and "or". So we extend the language with:

>           | Or [Prop] | And [Prop] deriving Show

The intention now is that And ... holds when all of the properties in the list hold and similarly for Or .... We can write an interpreter to test whether integers have the specified property:

> eval (Modulo k n) x = (x-k) `mod` n == 0
> eval (Or ps)      x = or $ map (\p -> eval p x) ps
> eval (And ps)     x = and $ map (\p -> eval p x) ps

(Note we limit ourselves to *finite* compositions of And and Or, otherwise eval wouldn't actually define a property due to non-termination.

There are lots of things we can say in our language. For example we can give the 'extreme' properties that are never true or always true:

> never = Or []
> always = And []

We can say that one number is divisible by another:

> divisibleBy k = Modulo 0 k

We can test it with expressions like:

*Main> eval (divisibleBy 3) 9
True
*Main> eval (divisibleBy 5) 11
False

We can also express non-divisibility. We say that n isn't divisble by k by saying that n is either 1, 2, ..., or k-1 modulo k:

> notDivisibleBy n =
>     let n' = abs n
>     in Or (map (\i -> Modulo i n') [1..(n'-1)])

(Disallowing n=0.)

*Main> eval (notDivisibleBy 3) 9
False

Eliminating And
It's not obvious at first sight, but there is a big redundancy in our language. There is no need for And. Consider And [Modulo k1 n1, Modulo k2 n2]. This asserts, for the number x, that x = s*n1+k1 and x = t*n2+k2. The Chinese remainder theorem tells us that either these have no solution, or that this pair of propositions is equivalent to one of the form x = k3 mod n3 for some k3 and n3. So every time we And a pair of propositions we can eliminate the And by using the theorem. Solving for k3 and n3 is straightforward. I use the extended Euclidean algorithm and the proof of the Chinese remainder theorem given at Cut the Knot.

> egcd n1 n2 | n2 == 0 = (1, 0, n1)
>            | n1 == 0 = (0, 1, n2)
>            | otherwise = (y, x-y*(n1 `quot` n2), g)
>            where (x, y, g) = egcd n2 (n1 `mod` n2)

> intersect (Modulo k1 n1) (Modulo k2 n2) =
>     let (s, _, g) = egcd n1 n2
>         (q, r) = (k2-k1) `quotRem` g
>     in if r == 0
>         then Modulo (q*s*n1+k1) (n1*n2 `quot` g)
>         else never

So now we can repeatedly use intersect pairwise on our properties to eliminate all uses of And. Here is some code to do so. Firstly, it's convenient to sometimes write any property as if it is a list of "subproperties", all Orred together:

> subproperties (Or ps) = ps
> subproperties p = [p]

Now we can go ahead and remove all of the Ands:

> removeAnd (Or ps) = Or (map removeAnd ps)

The property always can be rewritten as:

> removeAnd (And []) = Modulo 0 1

Remove And from the head of the list, remove it from the tail of the list, and then form all possible intersections of these two parts:

> removeAnd (And (p:ps)) = Or [q `intersect` q' |
>     q <- subproperties (removeAnd p),
>     q' <- subproperties (removeAnd (And ps))]
> removeAnd p = p

By induction, the return value from removeAnd can no longer contain an And. Note that the properties can grow in size considerably. Here is the proposition that x isn't divisble by 5 or 7 written out in full:

*Main> removeAnd (And [notDivisibleBy 5, notDivisibleBy 7])
Or [Modulo 1 35,Modulo 16 35,Modulo 31 35,Modulo 46 35,Modulo 61
35,Modulo 76 35,Modulo (-13) 35,Modulo 2 35,Modulo 17 35,Modulo 32
35,Modulo 47 35,Modulo 62 35,Modulo (-27) 35,Modulo (-12) 35,Modulo
3 35,Modulo 18 35,Modulo 33 35,Modulo 48 35,Modulo (-41) 35,Modulo
(-26) 35,Modulo (-11) 35,Modulo 4 35,Modulo 19 35,Modulo 34 35]

Now to the primes. Here's a standard way to make the list of primes in Haskell:

> isPrime primes n = foldr (\p r -> p*p > n || (rem n p /= 0 && r))
>                          True primes
> primes = 2 : filter (isPrime primes) [3..]

The Proof
Now we can give the proof this set is infinite. Suppose it were finite. Then we could form this property:

> prop = removeAnd $ And (map notDivisibleBy primes)

It contains no Ands, and so must simply be the Or of a bunch of Modulos. But each Modulo defines an infinite set, so prop must define an infinite set.

But prop is the property of not being divisible by any prime. So prop can only eval to True on -1 or 1, a finite set. Contradiction. Therefore primes is infinite.

We can look at approximations to prop like this:

> prop' n = removeAnd $ And (map notDivisibleBy (take n primes))

You can see that the proposition grows in size rapidly:

*Main> removeAnd (prop' 3)
Or [Modulo 1 30,Modulo (-83) 30,Modulo (-167) 30,Modulo (-251)
30,Modulo 71 30,Modulo (-13) 30,Modulo (-97) 30,Modulo (-181) 30]
*Main> removeAnd (prop' 4)
Or [Modulo 1 210,Modulo (-56159) 210,Modulo (-112319) 210,Modulo...]

Nonetheless, it would always be finite if there were only finitely many primes. As primes is infinite, you can think of the sequence prop' n as somehow trying to creep up on the set -1, 1, never quite getting there.

Unfortunately I have no time to explain why a topological proof should lead to one about a simple DSL beyond mentioning that there's a deeper story relating to the computability of eval for possibly infinite expressions of type Prop.

Addendum
I'll just say a little on the computability connection. Suppose we have a really dumb algorithm to test whether an integer x equals k mod n by doing a brute force search for s such that x=s*n+k. Suppose this is the only kind of test on x that we have available to us. The test will only terminate if it finds a solution. So with such an algorithm, testing for equality mod N is only semi-decidable. Now suppose we are allowed multi-threaded code. The infinitude of the primes implies that with our dumb tests, membership of -1,1 is also semi-decidable. So we can turn the problem of proving the infinitude of the primes into one about computability. You can see roughly how: we can "semi-test" Or ps by launching a process to test the first element of ps. Then launch a process to check the next, and so on. If any of these processes terminates, we have our answer. The argument presented above gives the details of how to construct a suitable Or ps.

An elementary way to approach Fourier transforms

2011-06-25T14:57:00.000-07:00

Why another introduction to the Fourier transform?
There are many elementary introductions to the discrete Fourier transform on the web. But this one is going to be a little different. I hope to motivate and demonstrate an application of the Fourier transform starting from no knowledge of the subject at all. But apart from this sentence, I'll make no mention of complex numbers or trigonometric functions. That may seem weird - the standard definitions seem to make it clear that these concepts are crucial parts of the definition. But I claim that in some ways these definitions miss the point. An analogy from software engineering is appropriate: most definitions of the Fourier transform are a bit like explaining an interface to an API by reference to the implementation. That might tell you how it works at a nuts-and-bolts level, but that can obscure what the API is actually for.

There's code all the way through so if anything I've said is ambiguous, that should help resolve it. I've chosen to write the following in 'literate Octave'. I can't say I love the language but it seems like the easiest way to express what I want here computationally.

Convolution
Suppose your camera has a hexagonal iris and it's out of focus. When you take a picture of a point light source (towards your lower left) you'll end up with a result like this:

The effect is known as bokeh and the wikipedia page has a nice example with an octagonal iris.

(If you want to use the Octave code in this article, save the image above as hexagon.png.)

Suppose we take a picture of an ordinary scene instead. We can think of every visible point in the scene as a point light source and so the out-of-focus photograph will be the sum of many hexagons, one for each point in the image, with each hexagon's brightness determined by the brightness of the point it comes from. You can also flip this around and think about the out-of-focus image as being a sum of lots of copies of the original image, one of each point in the hexagon. We can go ahead and code this up directly in octave.

Here's the picture that we'll apply bokeh to:

(Background: I used to work for a division of Kodak and that was the test image we frequently used. Save the image as marcie1.png.)

Let's read in our image and iris shape, converting the image to grayscale:

> I = double(imread('marcie1.png'))/256;
> I = (I(:,:,1)+I(:,:,2)+I(:,:,3))/3;
> J = double(imread('hexagon.png'))/256;


> h = size(I)(1);
> w = size(I)(2);

Now we'll apply the bokeh by looping over the entire iris, accumulating shifted copies of the original image. We're optimising a little bit. As many of the pixels in J are black we can skip over them. I'm using circshift so that we'll get wraparound at the edge. That turns out to be very convenient later.

> total = 0;
> K = zeros(h, w);
> for i = 1:h
>     i
>     for j = 1:w
>         if J(i, j) != 0
>             K = K + J(i, j)*circshift(I, [i, j]);
>             total = total + J(i, j);
>         endif
>     endfor
> endfor

We use total to scale the overall brightness back into a reasonable range:

> imshow(K/total)
> pause(5)

Make sure you understand what that code is doing because the rest of this article depends on it. The central line in the double loop is repeatedly adding copies of the original image, shifted by [i, j], and scaled by J(i, j).

There's just one problem with that code. It's incredibly slow. It's only tolerable because I know that most of J is zero so I could optimize it with a conditional. More general higher resolution images will leave you waiting for a long time.

The image we have computed above is known as the convolution of I and J. My goal is to show how we can use the Fourier transform to look at this in a very different way. As a side effect we will also get a much faster convolution algorithm - but the reason it runs faster is a story for another time. In this article just want to show what Fourier transforms are and why they're relevant at all.

Fourier transforms
Central to the code above is the fact that we're shifting and adding the same image over and over again. If we could avoid that we might find a way to speed the code up. So let me define some shorthand. I'll use R to represent the function that shifts an image right one pixel. In Octave that's given by circshift(I,[0,1]). I'll use U to mean a shift up by one pixel. Rather than define *the* Fourier transform I'm going to define a family of operations, all called Fourier transforms. (More precisely, these are two-dimensional discrete Fourier transforms.)

(1) A Fourier transform is a linear operation that converts an image to another one of the same dimensions. This means that if you apply it to sums and multiples of images you get sums and multiples of the Fourier transforms. In the language of Octave, if F is such an operation, then F(A+B) == F(A)+F(B) and F(s*A) == s*F(A), for s a scalar.

(2) A Fourier transform has this property: there is a pair of images, A and B, such that for any image I, F(U(I)) = AF(I) and F(R(I)) = BF(I). (AF(I) means the pixelwise product of A and F(I)). In Octave notation:

F(circshift(I, [1,0])) == A .* F(I)
F(circshift(I, [0,1])) == B .* F(I)

Fourier transforms convert shifts to multiplications. (If you only learn one thing from this article, this should be it.)

(3) Fourier transforms are invertible so there must be another linear transform F^-1 such that F^-1(F(I)) = I and F(F^-1(I)) = I.

Anything with these properties is an example of a Fourier transform. It's the second property that is crucially important. In jargon it's said to diagonalise translation.

From (2) it follows that we can compute the Fourier transform of an image shifted by [i, j] by multiplying the original Fourier transform by A.^i .* B.^j. (.^ is the pixelwise power function.)

If h is the image height, then shifting up h times should wrap us around to the beginning again. Similarly for w. So from (2) we know that A^h=1 and B^w=1. (That's shorthand for saying that each of the individual elements of A raised to the power of h give 1 and so on.)

It just so happens that Octave comes supplied with a suitable function that satisfies the three conditions I listed. It's called fft2. Let's find out what the corresponding images A and B are:

> r = rand(h, w);
> A = fft2(circshift(r, [1, 0])) ./ fft2(r);
> B = fft2(circshift(r, [0, 1])) ./ fft2(r);

(Note that ./ is the pixelwise division operator. If you're unlucky your random numbers will lead to a division by zero. Just try again.)

Let's try again for another random image:

> s = rand(h, w);
> A0 = fft2(circshift(s, [1, 0])) ./ fft2(s);
> B0 = fft2(circshift(s, [0, 1])) ./ fft2(s);

We can see that A is almost exactly A0 and B is almost B0:

> abs(A-A0)
> abs(B-B0)

So now we can go back to our original convolution algorithm. Let's define

> II = fft2(I);

And now we can use the first and third properties to compute the Fourier transform of the image with bokeh applied to it. We're applying properties (1) and (2) to the central line of the double loop above:

> KK = zeros(h, w);
> for i = 1:h
>     i
>     for j = 1:w
>         if J(i, j) != 0
>             KK = KK + J(i, j).*A.^i.*B.^j.*II;
>         endif
>     endfor
> endfor

I said that the definition of Fourier transform requires the existence of an inverse, and claimed that fft2 was a Fourier transform. So we must have an inverse. Octave conveniently provides us with the name ifft2:

> K = ifft2(KK/total);
> imshow(K)
> pause(5)

We've eliminated all that shifting, but at the end of the day this code is slower. But did you notice something curious about the innermost line of the double loop? It's always adding the same multiple of II to KK. We can completely factor it out. So we can rewrite the code as:

> JJ = zeros(h, w);
> for i = 1:h
>     i
>     for j = 1:w
>         if J(i, j) != 0
>             JJ = JJ + J(i, j).*A.^i.*B.^j;
>         endif
>     endfor
> endfor

We can leave the multiplication by II all the way to the end:

> KK = JJ .* II;
> K = ifft2(KK/total);
> imshow(K)
> pause(5)

This is pretty cool. We can precompute JJ. Any time we want to convolve an image with the hexagon we apply fft2 to the image, multiply by JJ and then apply ifft2. But there's something even better going on. Let's look more closely at that double loop above involving powers of the elements of A and B. Let's write out the function it computes in standard mathematical notation:

f(J) = Σ_i,jAⁱ B^j J_i,j

What is f(U(J))?

f(U(J)) = ΣAⁱ B^j J_i-1,j (dy definition of shifting up)
f(U(J)) = ΣAⁱ⁺¹ B^j J_i,j (by summing over i+1 instead of i and using wraparound)
f(U(J)) = A f(J)

Similarly,

f(R(J)) = B f(J)

In other words, f satisfies the second property of Fourier transforms. It obviously satisfies the first property. Our transform f looks like a Fourier transform. In fact, it is, and Octave's fft2 is defined this way. So now we have this procedure for convolving:

> II = fft2(I);
> JJ = fft2(J);
> KK = II .* JJ;
> K = ifft2(KK/total);

> imshow(K)
> pause(5)

We now have a fast convolution algorithm. Of course I've left out an important point: fft2 and ifft2 are fast. They don't use the obvious summation algorithm suggested by my definition of f. But that's an implementation detail. We're able to reason successfully about important properties of fft2 using the properties I listed above.

Conclusion
Let me recapitulate so you can see

1. I defined Fourier transforms
2. I showed how convolution could be rewritten using one
3. I told you that fft2 was a Fourier transform, giving an alternative algorithm for convolution
4. I showed that another part of the convolution algorithm looks like a Fourier transform. (I didn't prove it had an inverse.)
5. I told you that (by amazing coincidence!) this other part is also fft2.

I haven't shown you how to implement a Fourier transform. But I have shown you how you can reason about many of their properties. Enough to get much of the way to the convolution theorem.

Approaching Fourier transforms through the properties I listed is common in the more advanced mathematical literature. But for some reason, in the elementary literature people often choose to describe things in a more complicated way. This is true of many things in mathematics.

I hope the above serves as useful motivation when you come to check out a more standard and more complete introduction to Fourier transforms. In particular, now's a good time to try to understand why the usual definition of the discrete Fourier transform satisfies the properties above and so fill in the steps I missed out.

Note
You may need to set up Octave for graphics. On MacOSX I started the X server and used "export GNUTERM=x11" before running Octave.

Simulating visual artifacts with Fourier optics

2011-06-04T12:02:00.000-07:00

The problem

One of the fun things about working in the visual effects industry is the large number of disciplines involved, especially when working for a smaller company where everyone has to do everything. At one point many years ago we did some work on simulating camera artifacts including diffraction spikes. Although the wikipedia page is in the context of astronomy you can get diffraction artifacts with any kind of camera, and you'll get spikes whenever the iris has a polygonal shape. I hope to sketch why later.

I noticed an intriguing question about a diffraction pattern on Quora. (For those not following the link: this photo is a picture of a reflection of the photographer and camera in an LCD screen.) My immediate thought was "how could I have faked this in a render?".

My main interest was in the 'X' shape. And the first question is this: how do I know this is a diffraction effect? It might be that there are specular reflections off periodic structures in the LCD display orthogonal to the arms of the 'X'. This would explain everything with only geometric optics.

There are two big clues I know of that diffraction is at work: we have separation of white light into separate colours. This is typical of diffraction effects. We also see a discrete structure: rays directed along a discrete set of directions in addition to the arms of the 'X'. This is typical of diffraction from periodic microstructures, and again I hope to sketch why later.

So how can we explain the appearance of this photograph?

A simple model

If we consider diffraction effects we need to model light as waves. Let's start by working with just one frequency of light from the camera flash. Let's also suppose that the flash is far enough away from the monitor that the spherical wavefronts can be modelled as plane waves when they meet the monitor. Here's a picture:

The thick black vertical line is the monitor screen. The other vertical lines are the incoming plane waves. And the red arrow shows the direction of motion. Huygen's principle tells us that when the waves reflect off the monitor surface we can treat the entire monitor surface as a collection of points, each of which is emitting spherical waves equally in all directions. If you're unfamiliar with Huygen's principle it may seem wrong to you. When you shine a light at a mirror we get reflection along one particular direction, not in all directions. That's true, but we can view this fact as a consequence of interference between all of the little points on the surface emitting in all directions. Again, that's another thing that will emerge from the calculation below.

The question I now want to answer is this: how much light gets reflected in direction k? I'll assume we're observing things from far enough away that we can neglect the details of the geometry. We'll just compute the contribution in direction k (a unit vector (l, m, n)) from all of the points on our surface. We'll assume that because of variations on the surface, the intensity of the emission varies along the surface. We'll model the emission as a function I, so I(y, t) is the emission from position y at time t. We want to collect the light going in direction k through some light collector at time t. Call that J(k,t). Here's a picture:

Remember: there are waves going in all directions, but I've just drawn the waves going in direction k.

Fourier optics for beginners

Now comes the interesting bit: the light arriving at the collector is the integral of all the light emitted from the points on the monitor surface. But the time of flight from each point on the monitor surface is different. So when we perform our integral we need to build in a suitable delay. It's straightforward geometry to show that the time delay between waves arriving from y₁ and y₂ is my/c where c is the speed of light and m, as defined above, is the y-component of k. So, up to a constant of proportionality, and some absolute choice if time, we want this integral

J(k,t) = ∫ dy I(y,t-my/c)

We assumed that the monitor surface was being struck by plane waves. Assuming they were orthogonal to the surface and coherent this means that I is proportional to exp(iωt). The time dependence in J is just the sinusoid, so we drop that from further calculations. So we can write, ignoring another constant of proportionality,

J(k) = ∫ dy I(y)exp(-iωmy/c) = ∫ dy I(y)exp(-i2πmy/λ)

Where I used ω = 2πc/λ and λ is the wavelength of the light. In other words, J is the Fourier transform of I. More generally, if we modelled the z-axis we'd find that J is given by the 2D Fourier transform of I as a function of the 2D surface of the monitor. The actual power striking the sensor is given by the absolute value of the square of J.

Ordinary reflections

So now I can pop the first of my pending explanations off the stack. Suppose that the surface is completely uniform. Then I(y) = constant. The Fourier transform of a constant function is the dirac delta function. In other words - despite the fact that we're modelling every point on the surface as its own emitter, the resulting reflection from a perfectly smooth surface is a wave going straight back where it came from. (Exercise, modify the above treatment to show that waves striking a mirror at an angle obey the usual reflection law for mirrors.)

Periodic structure

Now suppose that the surface of our monitor has a periodic structure. The Fourier transform of a periodic function is concentrated at spikes corresponding to the fundamental frequency of the structure and its overtones. So we expect to see a grid-like structure in the result. I believe that's what we're seeing in the spray of quantised vectors coming from the center of the image and the fact that the arms of the 'X' look like discrete blobs of colour arranged in a line rather than the white spikes in the astronomical example I linked to above. That tells us we're probably seeing artifacts caused by the microstructure of the LCD pixels.

An edge

Consider the function I(y,z) = -1 for y<0 and 1 for y>0. In other words, the sgn function along the y-axis and constant along the z-axis. Up to multiplication by a constant, its Fourier transform is given by the dirac delta along the z-axis and 1/ω along the y-axis. On other words, it's a spike, starting at the origin, extending along the y-axis, but fading away as we approach infinity. This is the source of diffraction spikes. Any time we see such a spike it's likely there was a perpendicular edge somewhere in the function I. For telescopes it often comes from the struts supporting the secondary mirror. For cameras it comes from the polygonal shape of the iris. In this case, it looks like we must have pixels whose shape has edges perpendicular to the 'X'. I have to admit, when I came to this conclusion it sounded very implausible. But then someone posted this image from here:

There are those edges at the predicted angles.

Simulating the artifact

Now we have enough information to simulate the artifact. You can think of the remainder of this article as written in "literate octave".

Save the above close-up of the pixels in a file and read it in.

I = imread('lcd-pattern.jpg');

I crop out one quarter as it simply repeats. This also removes the text.

I1 = I(251:500,251:500,1);

Now clean up some stray (image) pixels and pick out just one (LCD) pixel.

I1(1:250,100:250) = 0;

The photo is blurry, probably has compression artifacts and seems quite noisy to me. So I threshold it at intensity 200 to get a hard shape. I then compute the discrete Fourier transform and scale the intensity to something nice.

I2 = 0.0001*abs(fft2(I1>200).^2);

The Fourier transform puts the origin in the corner so I shift it to the centre.

I2 = circshift(I2,[125,125]);

Now we need to do the same thing for red and green. We reuse the red image but it needs scaling in proportion to the wavelength. The numbers 1.2 and 1.4 are intended to be the ratio of the wavelengths to that of blue light. They're just ballpark figures: I chose them to make the image sizes convenient when rescaled. I also crop out a piece from the centre of the rescaled images so they line up nicely. I use imresize2 from here.

Ib = I2;

Ig = imresize2(Ir,1.2,'nearest')(26:275, 26:275);

Ir = imresize2(Ir,1.4,'nearest')(51:300, 51:300);

Now I assemble the final RGB image.

J = zeros(250, 250, 3);

J(:, :, 1) = Ir;

J(:, :, 2) = Ig;

J(:, :, 3) = Ib;

imwrite(J,'cross.jpg')

Here's the result:

Analysis

Not bad, but it's not perfect. That's unsurprising. Have a look at the thresholded image straight after thresholding. It's a very poor representation of the pixel shape. It would be better to redraw the shape correctly at a higher resolution (or even work analytically, not infeasible for simple shapes like these pixels). Stray (image) pixels can make a big difference in Fourier space. That explains the 'dirtiness' of my image. Nonetheless, it has the big X and it has the rows of dots near the middle. Qualitatively it's doing the right thing.

Note that the big X in the middle seems to have ghost reflections at the left and right. These are, I think, aliasing and purely artifacts of the digitisation.

It'd look slightly better with some multi-spectral rendering. I've treated the pixels as reflecting one wavelength but actually each one reflects a band.

There are also, almost certainly, lots more effects going in that picture. The close-up photograph is only revealing one part of the structure. I'm sure there is 3D structure to those pixels, not just a flat surface. I suspect that's where the horizontal white line is coming from. Because it's white it suggests an artifact due to something interacting equally with all visible wavelengths. Maybe a vertical edge between pixels.

But overall I think I've showed that Fourier optics is a step in the right direction for simulating this kind of effect.

Note

Everything I know about Fourier optics I learnt from my ex-colleague Oliver James.

All the errors are mine of course.

Fast forwarding lrand48()

2011-05-28T16:47:00.000-07:00

A break from abstract nonsense to answer a question I've seen asked online a number of times. It requires nothing more than elementary modular arithmetic and it ends in some exercises.

Given a pseudo-random number generator, say BSD Unix lrand48(), is there a quick way to jump forward a billion numbers in the sequence, say, without having to work through all of the intermediate numbers? The method is no secret, but I couldn't find explicit code online so I thought I'd put some here. Literate Haskell of course.

> {-# LANGUAGE ForeignFunctionInterface #-}
> {-# OPTIONS_GHC -fno-warn-missing-methods #-}

On MacOSX, if you type 'man lrand48', you'll see the function lrand48() returns a sequence of 31 bit non-negative integers defined using the sequence r_n+1 = ar_n+c mod m where

> a = 25214903917
> c = 11
> m = 2^48

The actual returned value is the floor of r_n/2¹⁷ and r₀ = 20017429951246.

We can compute the nth element in the sequence the hard way by importing lrand48 and looping n times:

> foreign import ccall "lrand48" lrand48 :: IO Int

> nthrand 1 = lrand48
> nthrand n = lrand48 >> nthrand (n-1)

But there is a better way. If we iterate twice we get that r_n+2 = a(ar_n+c)+c mod m = a²r_n+ac+c mod m. Note how two applications of the iteration give you back another iteration in the same form: a multiplication followed by an addition modulo m. We can abstract this a bit. Given two function f(x) = ax+c mod m and g(x) = a'x+c' mod m we get g(f(x)) = (a'*a)*x + a'*c+c' mod m. We can represent functions of this type using a simple Haskell type:

> data Affine = Affine { multiply :: Integer, add :: Integer } deriving (Show, Eq, Ord)

We can now write a function to compose these functions. I'm going to use the operator * to represent composition:

> instance Num Affine where
>    Affine a' c' * Affine a c = Affine (a'*a `mod` m) ((a'*c+c') `mod` m)

To skip forward n steps we just need to multiply n of these together, ie. raise Affine a c to the power of n using ^. We then need to apply this function to r₀:

> initial = Affine 0 20017429951246

> nthrand' n = (add $ Affine a c ^ n * initial) `div` (2^17)

Now try firing up ghci and comparing the outputs of nthrand 1000000 and nthrand' 1000000. Don't run nthrand more than once without resetting the seed, eg. by restarting ghci. (I know someone will post a reply below that it doesn't work...)

There are lots of papers on how to do this with other kinds of random number generator. My example is probably the easiest. The main application I can see is for jumping straight to that annoying regression test failure without going through all of the intermediates.

Exercises.
1. Read the corresponding man page for Linux. Port the above code to work there. Or any other OS you feel like. Or any other random number generator.
2. Can you split lrand48() into two? Ie. can you make two random generators that produce sequences s_i and t_i so that s₀, t₀, s₁, t₁, ... form the sequence given by lrand48().
3. I've neglected to mention some special sauce in the code above. Why does it actually run so fast? (Clue: why did I use Num?)

Perturbation confusion confusion

2011-04-30T10:34:00.000-07:00

Update: I'm making a bit of a turnabout here.

Firstly, I have to point out that at no point have I disagreed with S&P on any purely technical issue. We're all working with the same code and agree on what it does. The issue is a human one: is it possible to get wrong results by writing code that *looks* correct. I was sent this example:

d (\x -> (d (x*) 2)) 1

Read naively it produces a different result to what you might expect. We can see why it fails by looking at this:

d (\x -> (d (x*) 2)) (1::Integer)

It fails to typecheck! That "1" is actually of type D Integer. So the semantics of that code are entirely different to what you might expect if you read it like an ordinary mathematical expression and ignored the types of all the subexpressions.

So I agree that it is possible for a programmer to misread that code. I still don't consider the algorithm to have failed in this case. It is standard in Haskell that the appearance of an integer in code doesn't necessarily mean it's of Integer type, so that when we write Haskell code we always need to be aware of the types of all of our terms. When we write d (\x -> (d (x*) 2)) 1, we're asking for the wrong thing before the code has started executing. But I have been convinced that there are dangers here for people reading the code naively.

However I am now completely convinced the situation is more complex than this and I'll have to address it again.

Anyway, I suggest using AD regardless. It's an awesomely powerful technique. But don't capture variables in a lambda or function without wrapping them immediately them in a lift, and if you multiply nest, you need to multiply lift. Essentially lift is a signal to the compiler that you desire the variable to be held constant with respect to the derivative. It also makes the code more readable. It's analogous to using fmap the correct number of times in order to push a function down into nested lists.

So I'll leave the following in place because removing it would be bad form and some of it still stands.

Introduction
In a recent paper, Siskind and Pearlmutter warn how the use of differentiation operators in functional programming languages is "fraught with danger" and discuss a problem common to "all attempts to integrate a forward-mode AD operator into Haskell".

This is curious. I have had great success using automatic differentiation code both in Haskell and functional-style C++ and failed to notice the danger. Clearly I needed to take better notice of what I was doing.

So let's go ahead and implement AD and try to reproduce the problem they point out.

Automatic Differentiation

> data D a = D { real :: a, infinitesimal :: a } deriving (Eq, Show)

> instance Num a => Num (D a) where
>   fromInteger n = D (fromInteger n) 0
>   D a a'+D b b' = D (a+b) (a'+b')
>   D a a'*D b b' = D (a*b) (a*b'+a'*b)

We can now define a differentiation operator:

> d f x = infinitesimal (f (D x 1))

We can use d to differentiate a function like f(x) = x³+2x²+x+1 at 2 to get:

> example0 =  d (\x -> x^3+2*x^2+x+1) 2

Imagine you were confused enough by AD to write Siskind and Pearlmutter's example exactly as described in equation (2) of the paper:

example1 =  d (\x -> x*(d (\y -> x+y) 1)) 1

We don't get an incorrect result. Instead, we get this error message:

Occurs check: cannot construct the infinite type: a0 = D a0
Expected type: D (D a0)
  Actual type: D a0
In the first argument of `(+)', namely `x'
In the expression: x + y

The Haskell type checker identifies precisely where the problem is. x and y don't have the same types so they can't be added. Rather than immediately analyse how to fix this, let's try a completely different approach to calculus: symbolic differentiation. We'll define an expression type and write code to differentiate it

Symbolic Differentiation

> data E a = X | Const a | E a :+: E a | E a :*: E a deriving (Eq, Show)
> diff X = 1
> diff (Const _) = 0
> diff (a :+: b) = diff a + diff b
> diff (a :*: b) = a:*: diff b + diff a :*: b

We want to be able to evaluate these expressions:

> eval X x = x
> eval (Const a) x = a
> eval (a :+: b) x = eval a x + eval b x
> eval (a :*: b) x = eval a x * eval b x

We can make this easier to use by making E a an instance of Num:

> instance Num a => Num (E a) where
>   fromInteger n = Const (fromInteger n)
>   a + b = a :+: b
>   a * b = a :*: b

And now we can write an alternative to d:

> d' f x = eval (diff (f X)) x
> example1 =  d' (\x -> x^3+2*x^2+x+1) 2

We of course get the same result.

So let's try the example from the paper again:

example1 =  d' (\x -> x*(d' (\y -> x+y) 1)) 1

We don't get an incorrect result. Instead, we get this error message:

Occurs check: cannot construct the infinite type: t0 = E t0
Expected type: E (E t0)
  Actual type: E t0
In the first argument of `(+)', namely `x'
In the expression: x + y

An almost identical error message. So what's going on?

Look at the type signature to d':

d' :: Num t => (E a -> E t) -> t -> t

When we evaluate

d' (\x -> x*(d' (\y -> x+y) 1)) 1

the outer d' differentiates symbolically so internally it uses an expression type. This means that x is an expression, not a numerical value. But this means that the inner d' is being asked to evaluate its result at a value that is itself an expression, not a value. So internally it uses expressions of expressions. So x is of type E a and y is of type E (E a). It's no surprise we can't add them. Our bug is easily fixed. We define a lifting function:

> lift' x = Const x

and now we can correctly evaluate:

> example2 =  d' (\x -> x*(d' (\y -> lift' x+y) 1)) 1

> lift x = D x 0
> example3 =  d (\x -> x*(d (\y -> lift x+y) 1)) 1

Much the same discussion applies to the AD code. In that case, D a is the type of a's that have been infinitesimally perturbed. The type D (D a) is a value whose perturbation is itself perturbed. But Siskind and Pearlmutter say that we might be subject to perturbation confusion and could get 2. Curious. There is no way we can 'confuse' these distinct perturbations. Not no how. Not no way. They correspond to data of entirely different types. Far from being fraught with danger, the Haskell type system keeps us safe. Perturbation confusion in Haskell is about as likely as expression tree confusion.

Let's take a look at the analysis in section 2. They introduce a perturbation ε and show how we might have computed 2 instead of 1. But this derivation doesn't correspond to any computation we could possibly have made using the function d in Haskell. The only issue we have identified is that we have to write our code using correct types. This has nothing to do with automatic differentiation or perturbations. It applies equally well to the symbolic differentiation code. In fact, it has nothing to do with differentiation as it applies equally well to my code to compute the even and odd parts of a function.

We can press on with the paper. Section three tries to sell us a remedy - tagging. But we have no need for tagging. The Haskell compiler already deduced that the perturbation inside the inner d is of a different type to that in the outer d. The only explanation I can come up with for this section is that the authors have some experience with functional programming languages that are dynamically typed and are trying to apply that experience to Haskell.

Section 4 reiterates the point in the abstract that implementations of AD "fail to preserve referential transparency". I think we can safely ignore this claim. The AD code above isn't using any unsafe Haskell operations. It clearly is referentially transparent.

Confusion lifted
So now to section 6. They talk about the lift function and how correctly inserting it requires "sophisticated non-local analysis". Now everything becomes clear. Siskind and Pearlmutter don't consider this algorithm to be "automatic differentiation" unless they can leave out the lift operations. But it is not common practice to write Haskell code by deliberately writing a significant body of code that doesn't type check, and then expect to automatically insert the missing pieces. You simply write the code correctly to start with. In fact, when I first found myself inserting a lift in my own code I didn't even think of myself as having solved a problem - I just wrote the code that was analogous to code that Haskell programmers write every day. If I had forgotten the lift, the compiler would have set me straight anyway.

Conclusion
"Perturbation confusion" is a non-existent problem in Haskell AD. Siskind and Pearlmutter correctly point out that Haskell code using AD is not exactly analogous to mathematical notation for derivatives. This is unsurprising, differentiation is not computable (for deep reasons reasons that are similar to the reason that intermediate values are not computable). But functions like lift are an everyday part of typed functional programming. This is not a feature of AD, it applies equally as well to symbolic differentiation (and indeed many other parts of Haskell such as using monad transformers and nested functors). The algorithm is no less a form of AD because of the use of lift.

I have been advocating AD for many applications over recent years. It has been an uphill struggle. It is often incorrectly assumed that I am talking about numerical or symbolic differentiation. Programmers will often assert that what I am describing is completely impossible. But when they grasp what I am explaining the response often changes to "why didn't anyone tell me about this years ago?" It is disappointing to see two major contributors to computer science whose work I greatly respect scaring away potential users with FUD about a non-existent problem.

(And if you got here but skipped the update at the top you really need to go to the top as I now think Siskind and Pearlmutter have a good point! I'm leaving the above just for the sake of posterity.)

Image-based rendering and some ancient history

2011-04-09T09:37:00.000-07:00

Introduction
I've not said much about my work in visual effects in this blog. This is mainly because I try very carefully to avoid any kind of intellectual property conflict with my employer. But now I've left the world of visual effects to work at Google, I think I might throw in the occasional article. So here's something about my own work from many years ago. In particular, back in the year 2000 I and my colleagues George and Kim received an Academy Scientific and Technical Achievement Award for "the development of a system for image-based rendering allowing choreographed camera movements through computer graphic reconstructed sets." I rarely tell people what this work was for. But now I'll have a convenient web page to which I can direct anyone who asks. I'm hoping to aim this at people who know a little mathematics and geometry, but not necessarily anything about computer graphics.

Many years ago I worked at a company called Manex Visual Effects. One day the story of that company needs to be told. Here are some things that have been alleged about Manex that I neither confirm nor deny: its stock was traded in one of the largest ever pump-and-dump scams in the Far East, it extracted millions of dollars from Trenton, New Jersey for a fake scheme to create an East Coast "Hollywood", and it spent a couple of years making press releases about how it was doing work on the Matrix sequels despite the fact that there was nobody at the company who was actually doing work on the movies, including making press releases quoting me describing the ongoing work long after I had left. At one point I was apparently being trailed by a private detective who even came rowing a boat past the back of my house in order to get pictures of me. I haven't checked the fine print, but I think the contracts that prevented me from speaking about any of this expired long ago.

But back around 1998-2000, when Manex was doing real effects work, we developed a pipeline for a technique known as image based rendering. We became adept at taking large numbers of photographs of locations and then reproducing those locations as photorealistic renders. When I say photorealistic here I don't mean something that's merely supposed to look real, but actually looks fake. I mean renders that were indistinguishable from photography. That's commonplace now but back in 2000 it was a challenge, especially on a large scale.

The last ten seconds of this video clip from MI:2 should give some idea of what I'm talking about. The city around Tom Cruise as he jumps from the building is entirely CGI, but using real photography:

Texturing Cities
Although MI:2 isn't set in Sydney, that is the location that was used. A team returned with many hundreds, if not thousands of photographs of the central business district. We also managed to construct 3D models of these buildings by various means. We started by using Façade to construct the geometry, but ultimately we used a variety of methods including buying 3D models from, I think, the city itself. In later years I completely replaced the way we built geometry so we didn't need to use any other source.

The goal was to map the photography onto the models. Here's a picture illustrating three points being mapped onto a building. To render correctly we needed every visible point on every 3D object to be coloured (correct terminology: textured) using a pixel from one of the photographs.

The Matrix
We're ready for a theorem. Let P be a photograph of a scene S. Let (u, v) be ordinary rectilinear coordinates in P. Let (x, y, z) be 3D coordinates in the scene S. Define proj(u, v, w) = (u/w, v/w). Then there is a 3×4 matrix M such that for every point (x, y, z) visible in S, its colour is given by the colour of the point with coordinates proj(M(x,y,z,1)^T) in P.

(This assumes a pinhole projection model for the camera. In fact, real cameras have lens distortion. I wrote software to measure and fix this.)

The important point is that for each camera we just needed one matrix. We generated it by a very simple scheme: we had artists mark correspondences between points on our 3D models and points in the photography. I implemented a least squares solver to find the matrix that best fit each view. Pity the artists. They came fresh from college with well developed 3d modelling skills and I was responsible for reducing their careers to this simple point and click operation. But don't feel too bad. Most of these artists have since gone on to very successful careers in visual effects. They were a really great team.

But just mapping a single photograph is no good. The theorem tells us what to do for each point visible in a photograph, but how do we select which photograph to use? When producing the final render for production we had the luxury of being able to allow hours to render a single frame. We could make the decision on a pixel by pixel basis by performing a test for visibility from each camera, and then using heuristics to make the selection. (My colleague George had already proved this worked for the backgrounds to most of the bullet-time shots from The Matrix.)

But that was no good for the artists. They needed to see their work as they were doing it. If we could get the render time down to seconds it would be good. In fact, we got it down to a fraction of a second.

The Splitter
The first thing to note is that if we render a single triangle then the mapping defined by the theorem above is provided by modern graphics hardware. Back in 2000 it wasn't universally available, but it was in the SGI hardware we used. So we only needed to solve the camera selection problem fast. The solution was simple, we'd preprocess the geometry of the scene by splitting it into pieces (reducible to triangles), each piece corresponding a a part visible from one camera. Debevec et al. had already given demos of a technique for approximately partitioning a scene between cameras like this, but it didn't have the quality we wanted.

(By the way, my career in graphics started when I figured out an algorithm for eliminating the division from the proj() function and my new boss-to-be noticed I'd posted it online.)

There were three things needed:

1. For each camera, we needed to compute the subscene that was the intersection of the scene S with the volume viewed by the camera, its frustum.

2. For each subscene we needed to remove any parts occluded from view by any other part of the scene.

3. Ensuring that no two subscenes overlapped - ie. ensuring that no part of any subscene is visible in more than one camera.

It seemed like a messy computational geometry problem. But after a little thought it dawned on me that most of the geometry work, at least as measured by computation time, could be factored into one tiny geometry routine and the rest was logic. Here's a picture of the algorithm:

It takes a convex polygon and slices it using a plane (indicated in grey). Now everything else follows. For example, step 1 above can be achieved like this: the frustum associated with a camera is a 6-sided convex polyhedron (or a 4- or 5-sided polyhedron if you want it to extend all the way to infinity and/or all the way to the point of projection in the camera.) We can decide which points are in the frustum by slicing using the 6 planes and keeping the pieces that fall inside.

Step 2 works a little like this:

There is a 5-sided polygon lying in front of a box casting a shadow. The shadow volume is itself a 6-sided frustum (with a seventh side "at infinity", so to speak). So to remove the parts shadowed by a polygon we use the same slicing algorithm and collect up all the pieces that fall outside of the shadow volume. To remove all of the scene that is occluded from view by this camera we simply remove shadow volume corresponding to every single polygon in the scene. One of the virtues of image based rendering is that the detail in the photography makes up for the simplicity of the geometry, keeping the polygon count low. So the total number of operations might not seem too bad. Unfortunately, every time you slice through the scene you risk doubling the total number of polygons. It could have taken worse than exponential time in the number of polygons. But I took the risk, and surprisingly the time to preprocess was measured in minutes for typical scenes. The individual pieces of geometry were very intricate due to the fact that occlusion from any object could carve into a building. But like a jigsaw, the union of all the pieces gave something that looked just like a real city.

The last stage was ensuring that each part is visible from one camera. This was straightforward. As every polygon was processed with respect to every camera then I could associate to every polygon a list of cameras in which it was visible. At the end I could just sweep through and pick the best camera based on a heuristic. Typically we wanted to use the camera with the most undistorted detail.

Along the way I had to apply a few more heuristics to make things manageable. Many thin slivers would appear in the slicing. If they were small enough I threw them away. I'd also sweep through from time to time and fuse together neighbouring polygons that had been sliced but ended up still being visible in the same cameras, and whose union was still convex. That would reduce the total polygon count and speed up further processing.

It worked. The artists would mark their correspondences, run the optimiser to extract the per-camera matrices, kick off the 'splitter', have a coffee, and then return to a fully interactive view of Sydney, or wherever. They could now examine it from all angles and quickly ascertain what further work was needed, eg. if there were holes in one of the buildings. That allowed us to develop an efficient workflow and carry out the work on the scale needed to complete movie shots. We could also interactively visualise how other objects interacted with our cities. We could also use this tool to plan the photography and ensure we had all the coverage we needed.

In retrospect it all seems trivial. But at the time I don't think any other company could churn out entire city blocks the way we could. Today, Google do similar work on a much larger scale, and much more cleverly, with Street View.

And that was just one part of the technique that we marketed as "virtual cinematography" and which won George, Kim and I our award. But it's really important to remember that movie making is a big team effort. It took the work of many dozens of people to create our photorealistic cities, and the fact that I've chosen to write about my contribution doesn't mean that it was the only part that was important. It wouldn't have happened if George hadn't taught us the method originally, and if Kim hadn't believed in us and formed part of the company around our pipeline. And of course nothing wouldn't have happened if the artists didn't politely put up with our crudely engineered tools.

Ackowledgement
The example city image above was derived from a picture by Calvin Teo on wikipedia. I don't have the original photography we used.

Generalising Gödel's Theorem with Multiple Worlds. Part IV.

2011-04-02T15:04:00.000-07:00

Interpolating Between Propoitions
Suppose a, b and c are propositions. Then a∧b→b∨c. It seems that the a is irrelevant to what's happening on the right hand side. We could remove it and still have a true proposition: b→b∨c. It seems that the c on the right hand side is also irrelevant so that this is also true: a∧b→b. In fact, we can "refactor" the original proposition as a∧b→b→b∨c. (By Curry Howard it *is* a kind of refactoring.)

Let's try using the methods I described in Part I to demonstrate the validity of a∧b→b∨c:

I've done two things slightly differently. I brought the negation inside the implication at the start. This means we start with an implication with a clearly defined left- and right-hand side. It makes no difference to the outcome. But I've also coloured terms coming from the left-hand side in bLue and those on the right in oRange. Eventually we see a b and a ¬b meet each other resulting in the big red X telling us that the negation of a∧b→b∨c is invalid. But notice how the b and ¬b came from opposite sides of the implication. This tells us that the implication is valid because of the b's and that the appearance of a and c plays no role in establishing validity.

Now suppose a blue a met a blue ¬a. They would have both originated from the left hand side, telling us that the left hand side was a contradiction, regardless of what's on the right. So if the original implication we wished to establish were written as L→R then we'd have established that we could factor it as L→⊥→R. Similarly, if an orange a met an orange ¬a we'd have established the invalidity of ¬R meaning that we get L→⊤→R.

In all three cases we've managed to find an "interpolating" formula F such that L→F→R with the property that F only refers to letters that occur *both* in L and R. It may seem intuitively obvious that we can do this. Irrelevant hypotheses shouldn't play a role in establishing an implication. In fact, this is generally true of propositional calculus and is known as the Craig Interpolation Lemma. It also holds for Provability Logic. This is certainly not an obvious fact. The interpolation property fails for many logics.

I'm only going to roughly sketch how the proof of the interpolation lemma looks. Essentially it will be a proof by construction with the construction being the following code. We'll implement a bunch of tableau rules that construct the interpolation lemma. We've also seen what some of these rules look like in the discussion above.

Constructing Interpolations

I coloured propositions orange and blue above. I can't use colour in Haskell code so I'll instead label propositions L and R using the following type:

> data SignedProp a = L a | R a deriving (Eq, Ord, Show)

Think of the colour or sidedness of a proposition as its 'sign'. Sometimes we'll make a selection based on the sign:

> select (L _) l r = l
> select (R _) l r = r

> instance Functor SignedProp where
>     fmap f (L a) = L (f a)
>     fmap f (R a) = R (f a)

Sometimes we'll want to remove the sign:

> unsign (L a) = a
> unsign (R b) = b

Now you may see why I had to use Haskell's ViewPattern extension. I want to do the same pattern matching as before on these propositions even though they are a different type. As all the rules view the patterns through propType we can achieve this by making SignedProp a an instance of PropTypeable:

> instance PropTypeable a => PropTypeable (SignedProp a) where
>     propType (L a)    = fmap L (propType a)
>     propType (R a)    = fmap R (propType a)
>     neg               = fmap neg
>     isF               = isF . unsign
>     positiveComponent = positiveComponent . unsign
>     negative          = negative . unsign

With this framework in place the code to interpolate is quite small. Instead of returning a Bool or a diagram, these rules return Just the interpolating proposition when possible, otherwise a Nothing. But the code is also going to do something slightly more general. If we close a potential world containing propositions l₁,l₂,...,r₁,r₂,... we have shown that l₁∧l₂∧...∧r₁∧r₂∧... isn't valid. Ie. that l₁∧l₂∧...→¬r₁∨¬r₂∨... *is* valid. If the l_i are L-propositions, and the r_i are R-propositions, then the rules I define below will find an interpolating formula for l₁∧l₂∧...→¬r₁∨¬r₂∨.... For the original case of p→q we prime it with p on the left and ¬q on the right:

> interp (p :-> q) = let t = runTableau interpRules [L p, R (Neg q)]
>                    in simplify <$> t

And now we can give the rules:

> interpRules = TableauRules {

The first three branches of this case are the three rules discussed at the beginning. The fourth should be fairly obvious:

>     foundContradiction = \a -> Just $ case a of
>        (L _, L _) -> F
>        (R _, R _) -> T
>        (L n, R _) -> n
>        (R n, L _) -> Neg n,

The algorithm needs to know when a world has closed. The original validity rules returned a boolean. These rules return a Maybe, and we know a world closed if the algorithm succeeded in returning an interpolating proposition:

>     closes = isJust,

These are the trivial cases of finding a blue or orange ⊥:

>     foundF = \a -> Just $ select a F T,

Some rules that really just serve as glue:

>     open     = \_ -> Nothing,
>     doubleNegation = \_ t -> t,
>     conjRule = \_ _ t -> t,

Now suppose we want to find an interpolating formula for l₁∨l₂→r. If we set up our tableau we get:

(Note I'm using l and r as metavariables here, so l₁, l₂ and r represent propositions made up of (possibly many) ordinary single letter variables.) If we complete the two sides of the divide using our interpolation algorithm recursively we'll find propositions f_i such that l_i→f_i→r. Hence we find l₁∨l₂→f₁∨f₂→r. Clearly the middle proposition only contains letters that appear on both sides of the original implication. A similar analysis allows us to find an interpolating proposition for l→r₁∧r₂. That gives us this rule:

>     disjRule = \p _ _ tl tr -> select p (:\/) (:/\) <$> tl <*> tr,

Those are the rules for classical propositional calculus.

Now comes the tricky bit. If we draw a big red X in a subworld it allows us to back out and deduce an interpolating proposition in the parent world. The rule is simple but I'll leave the proof sketch as an appendix:

>     processWorld  = \p t -> select p Dia Box <$> t,

If any subworld of a world is closed, so is the parent world, so we can ignore all but the first closed subworld:

>     combineWorlds = mplus,

>     tableau = \_ t -> t
> }

You can reproduce the wikipedia example with this interp ((neg (p /\q) --> neg r /\q) --> (t --> p) \/ (t --> neg r)).

Definability

In my previous post I talked about how in provability logic we can define propositions implicitly. Now we're in a position to do the construction. Firstly we need to say precisely what we mean by a definition inside the language of provability logic, and then we need to say what ensures that definitions make sense.

An implicit definition of a proposition in provability logic is a function Prop -> Prop that doesn't analyze its argument. Informally, it defines a proposition p if f p is valid. Some candidates might be:

> def1 p = T --> p
> def2 p = p --> T

The first, def1, seems fine. It seems to uniquely single out p == T because def1 T is valid. But def1 (neg F) is also valid. So we can't uniquely pin down propositions, but only up to some sort of equivalence. In fact, we can use <-> as our equivalence relation.

But the second attempted definition is useless. Any proposition satisfies it. So we only consider a definition d to be valid if any two propositions satisfying it are equivalent. So we can say that d defines a proposition h if d p --> (p <-> h) is valid.

So here's the beginning of a function that attempts to find a proposition satisfying an implicit definition:

> beth d = let p = Letter "__p"
>              q = Letter "__q"
>          in if not $ valid $ (d p /\ d q) --> (p <-> q)
>             then error $ show (d p) ++ " doesn't satisfy precondition"

You can see how I've encoded our precondition for a definition to be good. I also used __p and __q to make sure we didn't clash with any letters in the definition. It's called beth because the theorem that ensures we can write the else clause is known as the Beth Definability Theorem. The last line is astonishingly short:

>             else interp (d p /\ p --> (d q --> q))

Suppose we have proved that d(p)∧d(q)→(p↔q). It immediately follows that d(p)∧p→(d(q)→q). (It's essentially a bit of currying.) This is a candidate for Craig interpolation as the left hand side has no q and the right hand side has no p. So we can make a sentence h so that d(p)∧p→h and h→(d(q)→q)). The lettes p and q are just letters. If a proposition is true for q, it's also true for p. So with a little rearrangement we get that d(p)→(p↔h). Or in English, if p satisfies our definition, h is equivalent to it. So we've constructed an h that does what we want.

Fixed Points
Now it's a small step to get a fixed point. We just make a definition of fixed point and apply beth. This looks like it might work:

> isFixedPoint' f p = p <-> f p

Unfortunately it's possible for p and p' to be inequivalent and yet both satisfy isFixedPoint' f. We have to use a "stronger" definition. I'll leave the proof to Boolos's book, but what we'll do is assert not just that p is a fixed point, but also that it is provably so. So we use:

> box' a = a /\ Box a
> isFixedPoint f p = box' (p <-> f p)

In order for the uniqueness to work, every occurence of the argument of f must be inside a Box or Dia.

And that's it. We can now churn out fixed points to our heart's content.

Gödel's Second Incompleteness Theorem again
Let's find a proposition that asserts its own unprovability:

> godel = fixedpoint $ \p -> neg (Box p)

We get Dia T, which is the claim that arithmetic is consistent. So godel shows that if arithmetic is consistent, then it can't prove its consistency.

Here's a large number that I use for regression testing that I lifted from Boolos's book:

> fpexamples = [
>         (\p -> Neg (Box p), Dia T),
>         (\p -> Box p, T),
>         (\p -> Box (Neg p), Box F),
>         (\p -> Neg (Box (Neg p)), F),
>         (\p -> Neg (Box (Box p)), Dia (Dia T)),
>         (\p -> Box p :-> Box (Neg p), Dia (Dia T) \/ Box F),
>         (\p -> Box (Neg p :-> Box F) :-> Box (p :-> Box F),
>                Dia (Dia (Neg (Box F) /\ Neg (Box F)) /\
>                Neg (Box F)) \/ Box (Box F)),
>         (\p -> Box p :-> q, Dia (Neg q) \/ q /\ q),
>         (\p -> Box (p :-> q), Box q),
>         (\p -> Box p /\ q, Box q /\ q),
>         (\p -> Box (p /\ q), Box q),
>         (\p -> q \/ Box p, T),
>         (\p -> Neg (Box (q :-> p)), Dia q),
>         (\p -> Box (p :-> q) :-> Box (Neg p),
>                Dia (Box F /\ Neg q) \/ Box F),
>         (\p -> q /\ (Box (p :-> q) :-> Box (Neg p)), q /\ Box (Neg q)),
>         (\p -> Dia p :-> (q /\ Neg (Box (p :-> q))),
>                Box F /\ Box F \/ q /\ Dia ((Box F /\ Box F) /\ Neg q)),
>         (\p -> Box (Box (p /\ q) /\ Box (p /\ r)),
>                Box (Box q /\ Box r))
>     ]

Note how in every case, p is inside a Box or Dia. It's pretty mind-bending to try to think about what all of these propositions could possibly mean.

We can easily write to test to see whether two propositions are equivalent (in the sense that they imply each other):

> equiv p q = valid (p <-> q)

> regress2 = do

>     print $ and $ map (\(f, x) -> fromJust (fixedpoint f) `equiv` x) fpexamples

I am disconcerted that one of the easy looking examples in Boolos fails my tests. Given that the difficult cases agree with my code I think it is likely an error by Boolos, though that's a scary claim to make against one of the best known logicians in the world.

The company Theory Mine has been selling theorems. But they don't seem too interesting. Far better to generate theorems that generalise Gödel's work and tell you about the very nature of provability. Just send me $19.99 and I'll send you a certificate.

> main = do
>           regress1
>           regress2

Notes
Here is some Prolog code by Melvyn Fitting to do much the same thing. Note that code relies crucially on backtracking whereas my code explicitly searches through subworlds.

The function runTableau is a kind of generalised fold. Such things have associated induction principles. In this case it's a generalisation of the "Unifying Principle" defined by Smullyan in First-Order Logic. I guess it is an ordinary generalised fold over a tree structure representing fully expanded tableaux, but we never explicitly build such a structure.

Oh...and apologies (if anyone has come this far) for how long this took. Last year I threw together code to find fixed points in a couple of hours thinking "that was easy, I can blog about it". But I hadn't realised how many hours of work it would take to explain what I had done. And even if nobody came this far, the rubber ducking improved my own understanding greatly.

Appendix
Suppose we successfully showed that this closed:

Then the generalised Craig Interpolation Lemma says that for some suitable f, both of the following close:

The way we'll prove this is to follow what we must have done to get the first potential world to close, assume inductively that the interpolation lemma works for all of the potential subworlds that we entered.

So let's start by supposing that we used ◊l₂ to open up a subworld, giving us:

If the subworld closed then recursively using interpolation we know that for some appropriate choice of g, these closed too:

Now consider trying to close these two worlds:

Using Worlds 1 and 2 above we find they close immediately:

And that now tells us that when we recurse back up, we can use f=◊g.

I started with world containing just 6 propositions so I seem only to have proved this for that case. But in fact, any combination of propositions just ends up with an argument that is substantially the same. The only thing that might change is that if we use ◊r₂ to try to open up a subworld we find that f=◻g for some other g. And this is all summarised compactly by the rule processWorld above.

Build Yourself a Bluetooth Controlled Six-Legged Robot

2011-02-27T14:42:00.000-08:00

Introduction

I don't know why you can't easily find toy robots that can be controlled by Bluetooth. But that lack can be remedied by a few hours work. So you can see exactly what I mean, here is the robot I'm going to describe being controlled by me via a program running on a Mac:

As it's controlled from a much more powerful computer you can use sophisticated algorithms to direct it. And if you don't want to do that, it can easily be programmed to run autonomously. Either way, you can make use of a pair of infrared proximity sensors as input to your algorithms.

The first thing I must do is give credit where it's due. This is a modification of someone else's design. In fact, it's a sample project at Pololu Electronics and Robotics modified by addition of an off-the-shelf Bluetooth board.

You'll need:

Everything required to build the Pololu sample project. That includes Windows. I was developing on a Mac so I made use of VMWare. I can't vouch for any other Windows emulation because you'll need a USB port and I've experienced troubles with emulation and USB in the past. You'll only need Window for initial configuration of the robot. After that you can use the development platform of your choice. The list of parts is here. Order the partial kit for the motor controller. You'll need the header pins it comes with.
One BlueSMiRF Gold Bluetooth modem from Sparkfun.

The first step is to build the Pololu project exactly as described. At this point, you'll have a complete robot that can run autonomously. You can program it from Windows using Pololu's scripting language described in the manual.

Note that the battery used is this. Its connector is different to that described in the project. You simply need to solder a pair of pins rather than a special type of connector:

By the way, there is a minor flaw in Pololu's design. The sensor and servos are glued to the battery using hot glue. When the battery is charged it heats. You can guess the rest. So once everything works, you may want to switch to a different glue. I used epoxy.

The strands of wire from the servos are very fine. At one point I think some came astray causing a short-circuit. I was surprised the controller board still worked after I saw smoke rise up from it. So after that I used some dabs of hot glue to act as insulation around my soldering.

I found that I couldn't use the motor controller via an external USB hub, I always had to use a cable directly from my Mac.

In the process of making the robot you cut the connectors off the servos. Keep two of these with 2 or 3 inches of wire attached. We'll reuse them later.

Setting up Bluetooth

In its default mode, the robot is scriptable and controllable via the USB port. However, two pins on the board are attached to a serial port (running at 5V, not the usual RS232 voltage). These can be connected directly to a pair of pins on the BlueSMiRF board.

First we'll get the modem working independently of the robot. So go ahead and solder 6 header pins to the board like so:

Now we can connect it directly to the robot's battery using its connector to make the following circuit:

Conveniently, VCC and GND are next to each other, just like in the battery connector. The little red LED marked Stat should start flashing. This means it's ready to go. Take care not to reverse the polarity (this isn't Star Trek) or make an off by one pin error.

Now go ahead and pair it with your main computer. On the Mac it's the usual process. It'll probably appear as a device called Firefly-something. There have been various revisions of documentation for BlueSMiRFs. I guessed, correctly, that the correct passkey was 1234.

Once it's connected you should get a new serial port. Mine appeared as /dev/tty.FireFly-DB29-SPP. SPP is Serial Port Profile. Now you need a terminal application that can talk to a serial port. On Unix machines you can use screen followed by the path to the dev file. You'll find help on the web for other OSes. On Windows I expect you'll get a new COM port. It didn't seem to matter what baud rate I chose at this stage. I guess that SPP presents an interface that looks like a serial port but that the connection rate doesn't really mean anything. But I may be wrong and you may need to experiment and/or read documentation. It may take several attempt to successfully connect. (It likes to prove it has the power to say no if it wants.) When you do, you'll see the Conn LED light up green.

If you connect within one minute of powering up the BlueSMiRF you can configure it by entering its command mode. One revision of the documentation says you should type +++. My version used $$$. You won't see your keys echoed on the screen but it should reply CMD. Now type d and return to get a summary of its status. Now type SU,19200 (and return) to set the baud rate on the serial port side. It'll go much faster than this but it's worth being cautious to start with. It acknowledges this but it doesn't actually change speed until you restart the device (eg. by power cycling it). Quit command mode by typing ---. At this point, anything you type in the terminal will come out at 19200 baud on the TX pin. Quit screen by hitting ^A^X.

I'm not sure how you can easily test this short of completing the robot. I tested it using a program on a MSP430 Launchpad. I wish I had a fancy diagnostic device to read serial transmissions. I ought to build one.

By the way, talking to the BlueSMiRF as if it's a serial port is a form of backward compatibility. It should be possible to talk to it directly via the RFCOMM API. Among other things, this would allow better diagnostics.

Configuring the Robot

Now you need to configure the robot to use 19200 baud on its TX/RX pins. Connect via USB to the Windows Control Center application (like you've done before if you built the original design) and set it to use the UART at 19200 baud. Make sure CRC is off. Like so:

Once that is done you don't need Windows again unless you change the baud rate.

Connecting Modem to Robot

Now you need to connect the robot to the modem. You have 6 pins on the modem board so you need a 6 pin connector. I glued together the two three pin connectors I mentioned above to make a connector like this:

The two outer wires can be soldered together because we're not using flow control, and so the BlueSMiRF is clear to send (CTS) whenever it is ready to send (RTS), The ends of the other four wires can be soldered directly to the motor controller board to make this circuit:

Now comes the iffy bit. I hot glued the BlueSMiRF board directly to the side of the servos like this:

That's iffy because it places electric motors directly next to a device using RF. Motors emit lots of RF noise. But they can't be that bad, after all servos have been used for years on radio controlled models, and Bluetooth probably has some error correction. It seemed to work for me. (But see note below.)

Talking to the Robot

Now we need to talk to the robot. Just about every development platform has a way of talking to serial ports. I used Haskell with the serialport package. As I had fixed the baud rate, I used the Pololu compact protocol. As the focus in this article is on hardware, I'll just give the code at the end. Once it's running you should be able to control the robot using

f - take a step forward
b - take a step back
l - turn left
r - turn right
e - read the state of the two proximity sensors

Now you're free to code up whatever you want. I experimented with localization using a method I learnt from Eric Kidd.

Notes

I occasionally had dropped bytes. I haven't yet worked out a way to 100% reliably recover when this happens and I'm not sure if it's due to the RF noise I mentioned above. You can probably mess with the timeout parameter in the Pololu Control Center to ensure the robot returns to a known state if communications cease for a bit.

The Apple SPP driver is pretty crappy. It generally works (apart from the annoying need to make multiple attempts to connect), but if you interrupt your code at the wrong time with ^C (say) you can find it locks up so hard that you end up with an unkillable zombie process. The OS didn't even seem able to kill it on shutdown (maybe it would have eventually timed out) and I had to power down the Mac if I wanted to use the robot again. It only happened a couple of times in two weeks of intensive use, but it's annoying.

And if the robot stops working - remember these batteries are pretty small and are doing a lot, so you may just need a recharge.

Code

The code uses two threads, a server that talks to the serial port and a client for users. This means that independent threads can control the legs and eyes, say, without the serial port transactions becoming interleaved. This isn't very polished but should be good enough to start experimentation. I've only tested it under Mac OS X but it looks platform independent to me.

> module Robot where

> import Prelude hiding (Left, Right)
> import Control.Concurrent
> import Control.Monad.Trans
> import Control.Monad.Trans.Maybe
> import Control.Monad.Trans.Reader
> import Control.Exception
> import System.Hardware.Serialport

Server side using Pololu compact protocol.

> sendByte :: Int -> SerialPort -> IO ()
> sendByte b s = sendChar s (toEnum b)

> getByte :: SerialPort -> IO (Maybe Int)
> getByte s = do
> b <- recvChar s
> return $ fmap fromEnum b

> setTarget :: Int -> Int -> SerialPort -> IO ()
> setTarget channel target port = do
> sendByte 0x84 port
> sendByte channel port
> let (b1, b0) = target `divMod` 0x80
> sendByte b0 port
> sendByte b1 port

> getPosition channel port = runMaybeT $ do
> liftIO $ sendByte 0x90 port
> liftIO $ sendByte channel port
> lo <- MaybeT $ liftIO $ getByte port
> hi <- MaybeT $ liftIO $ getByte port
> return $ lo + 0x100*hi

> data SerialCommand = SetTarget Int Int
> | GetPosition Int (Maybe Int -> IO ())
> | End

> serialExec command port =
> case command of
> SetTarget channel target -> setTarget channel target port
> GetPosition channel continuation ->
> getPosition channel port >>= continuation
> otherwise -> return ()

> serialThread commandChannel port = do
> command <- liftIO $ readChan commandChannel
> serialExec command port
> case command of
> End -> do
> closeSerial port
> return ()
> otherwise -> serialThread commandChannel port

Client side

> data Servo = Servo { channel :: Int, loLimit :: Int, hiLimit :: Int }

> data Direction = Left | Right | Forward | Back

> start tty = try (openSerial tty defaultSerialSettings
> { baudRate = B19200 }) >>=
> either
> (\ex -> do
> print (ex :: IOException)
> threadDelay 250000
> start tty)
> (\port -> do
> commandChannel <- newChan
> forkIO $ serialThread commandChannel port
> return commandChannel)

> type RobotM a = ReaderT (Chan SerialCommand) IO a

> writeChan' = flip writeChan

> end = ReaderT $ writeChan' End

> setServo limit servo = ReaderT $ writeChan' $
> SetTarget (channel servo) (limit servo)

> setLo = setServo loLimit
> setHi = setServo hiLimit

> readEye channel cont = ReaderT $ writeChan' $ GetPosition channel cont

> delay = 100*1000

> raise Left = setLo midLegs
> raise Right = setHi midLegs
> forward Left = setLo leftLegs
> forward Right = setHi rightLegs
> back Left = setHi leftLegs
> back Right = setLo rightLegs

> opposite Left = Right
> opposite Right = Left
> opposite Forward = Back
> opposite Back = Forward

> move Forward = forward
> move Back = back

> halfCycle side direction0 direction1 = do
> raise side
> liftIO $ threadDelay delay

> move direction0 Left
> move direction1 Right
> liftIO $ threadDelay delay

> motion direction0 direction1 = do
> halfCycle Left direction0 direction1
> halfCycle Right (opposite direction0) (opposite direction1)

> walkCycle = motion Forward Back
> reverseCycle = motion Back Forward
> turn Right = motion Forward Forward
> turn Left = motion Back Back

> withRobot tty cmds = do
> commandChannel <- start tty
> flip runReaderT commandChannel (cmds >> end)

> eyes = do
> readEye 4 $ print . (("Left eye = " ++) . show . fmap (<512))
> readEye 3 $ print . (("Right eye = " ++) . show . fmap (<512))

> commandLoop = do
> key <- liftIO $ getChar
> case key of
> 'f' -> walkCycle
> 'b' -> reverseCycle
> 'l' -> turn Left
> 'r' -> turn Right
> 'e' -> eyes
> otherwise -> return ()
> if key=='q'
> then return ()
> else commandLoop

> testRobot = do
> withRobot tty $ do
> liftIO $ print "Ready"
> commandLoop

Here's some stuff for you to configure. In particular, set the limits of the servos so that the legs don't whack into the body of the robot, stripping the gears. These numbers will be a little different depending on exactly how you built the robot. You may want to start with the upper and lower limits closer to 6000, the middle of the range of the servo motion.

> tty = "/dev/tty.FireFly-DB29-SPP"

> rightLegs = Servo 0 5000 6500
> midLegs = Servo 1 5000 6500
> leftLegs = Servo 2 5000 6500

One Last Thing
Double check everything I say above. I could have easily made a mistake. In particular, make sure that I haven't inadvertently introduced short circuits by comparing my diagrams against the online documentation for the parts.

Generalising Gödel's Theorem with Multiple Worlds. Part III.

2011-02-25T15:55:00.000-08:00

Fixed Points
A number of branches of mathematics have some sort of implicit function theorem. These guarantee that we can define functions implicitly, rather than explicitly. For example, we can define the function f from the positive reals to the positive reals by the relation f(x)^2 = x. In this case we can write an explicit formula, f(x) = x^1/2, but the implicit function theorem gives quite general conditions when the implicit equation defines a function uniquely, even when it is very hard to write an explicit formula. We have something similar in the theory of datatypes. We define the list type implicitly in terms of lists: L(X) = 1+X L(X). There is a unique smallest type (and unique largest type) that satisfies this equation and we use such equations in Haskell programs to uniquely specify our types. (Haskell uses the largest type.)

Provability logic also has a kind of implicit function theorem. Consider the following function:

> f0 p = Neg (Box p)

Suppose the equation p <-> f0 p had a solution. Then we would have a proposition that asserts its own unprovability. In fact, we know such a proposition: the Godel sentence from Godel's first incompleteness theorem. However, the Godel sentence uses a clever Quining technique to make a proposition that makes assertions about its own Godel number. We can't express such a thing in the language of provability logic. However, consider Godel's second incompleteness theorem. This tells us that if PA can prove its consistency then it is in fact inconsistent. Or to turn it around, if PA is consistent, then PA can't prove its consistency. Asserting that PA is consistent is the same as ◊⊤ (because asserting that ⊤ is consistent with the rest of PA is the same as saying the rest of PA is consistent). So Godel's second incompleteness theorem can be written in provability logic as ¬◻(◊⊤)→◊⊤. Using the methods of part one we can show that ◊⊤→¬◻(◊⊤) is also valid. So Dia T is the solution to p <-> f0 p:

> test0 = valid $ let p = Dia T in p <-> f0 p

Considered as a data structure, notice how f0 places its argument inside a Box. Consider just the functions Prop -> Prop that (1) simply copy p into various places in its return value (ie. that don't analyse p in any way) and that (2) put all copies of its arguments somewhere inside a Box (possibly in a deeply nested way). We'll call these 'modalising' functions. Then there is an amazing result due to Solovay:

The Fixed Point Theorem
If f is modalizing then it has a unique fixed point in provability logic. (Unique in the sense that if p and q are fixed points for f, then p <-> q is valid.

Our goal in the next post in this series will be to implement an algorithm to find these fixed points.

Think about what this means. We can make up any old function that has this property and find a corresponding Godel-like theorem. For example, pick

> f1 p = Box p --> Box (Neg p)

This has a fixed point p = Box (Box F) -> Box F which we can test with

> test1 = valid $ let p = Box (Box F) --> Box F in p <-> f1 p

Or in English (quoting verbatim from Boolos) "a sentence is equivalent to the assertion that it is disprovable-if-provable if and only if it is equivalent to the assertion that arithmetic is inconsistent if the inconsistency of arithmetic is provable". (That's a good argument for using provability logic instead of English!) These are the generalisations of Godel's theorem in my title.

But there are many other important consequences. The solution to this fixed point equation doesn't involve any self-reference (because we can't do self-reference in provability logic as we have no way to talk about Godel numbers). Self-reference is less essential than you might think - we can solve these equations without using it.

At first the fixed point theorem seems like a positive thing: we can crank out as many of these theorems as we like. But there's also a flip side: it says that once we've proved Lob's theorem, there aren't any more techniques we need in order to construct fixed points. So there aren't any funky new variations on what Godel did, waiting to be discovered, that will give us a bunch of new fixed points. We have them all already.

In order to find these fixed points I'm going to use the runTableau function we wrote last time. But that's for the next article in the series. In this post I want to put it to a simpler use:

Automatically Drawing Tableau Diagrams
First we'll need a very simple library of functions for drawing ASCII art diagrams. We'll represent a drawing simply as a list of ASCII strings, one for each row. Think of diagrams as forming a box. The width of a box is the length of its longest row:

> width box = foldr max 0 (map length box)

Sometimes we'll want to pad the length of the rows so that they all have the same length:

> pad len b a = a ++ replicate (len - length a) b

The aside function allows us to 'typeset' one box next to another:

> aside a b = let w = width a
>                 h = max (length a) (length b)
>                 a' = pad h [] a
>                 b' = pad h [] b
>                 a'' = map (pad w ' ') a'
>             in zipWith (++) a'' b'

We can draw nice frames around our boxes:

> frame a = let w = width a
>               h = length a
>               strut = replicate h "|"
>               rule = "+" ++ replicate w '-' ++ "+"
>           in [rule] ++ strut `aside` a `aside` strut ++ [rule]

And the rule for disjunctions requires us to draw side by side boxes with a vertical line between them:

> alt a b = let h = max (length a) (length b)
>               strut = replicate h " | "
>           in a `aside` strut `aside` b

Generating diagrams is now a matter of generating a diagram for each of the 'hooks' in our algorithm:

> diagramRules = TableauRules {

We'll modify the rules so that instead of simply returning a Bool to indicate closure we return a pair. The first element is a list of strings representing the rows of the ASCII art, but the second element plays the same role as before:

>     closes             = snd,

These are the functions that handle discovery of a contradiction. They draw a "X" in a circle to indicate this:

>     foundF             = \_ -> (["(X)"], True),
>     foundContradiction = \_ -> (["(X)"], True),

>     open               = \_ -> ([], False),

We deal with conjunctions simply by listing the two subpropositions that went into them:

>     conjRule           = \a b -> first ([show a, show b] ++),

For disjunctions we draw the diagrams:

>     disjRule           = \p a b (tl, lb) (tr, rb) ->
>                             (([show a] ++ tl) `alt` ([show b] ++ tr), lb && rb),

>     doubleNegation     = \q -> first ([show q] ++),
>     combineWorlds      = \(t0, b0) (t1, b1) -> (t0 ++ t1, b0 || b1),
>     processWorld       = \_ t -> t,

When we create a box for a new world we 'import' some propositions from the parent world. We mark these with an asterisk:

>     tableau            = \br (t, b) -> (frame (map (("* " ++) . show) br ++ map ("  "++) t), b)
> }

> diagram p = do
>     let (t, b) = runTableau diagramRules [Neg p]
>     mapM_ putStrLn t
>     print b

Now we can prove that we really do have a fixed point of f1:

> diagram1 = diagram $ let p = Box (Box F) --> Box F in p <-> f1 p

You may want to pick a really small font and stretch your terminal wide before you run that:

See how every box ends in an (X), just as we want.

note
To run the code above, just append this article to the previous one to make a single literate Haskell program.

Quine Central

2011-01-30T18:01:00.000-08:00

Quines for Everyone
This is an interruption to the sequence on provability logic (though it's not entirely unrelated).

The code below spits out a Haskell program that prints out a Perl program that prints out a Python program that prints out a Ruby program that prints out a C program that prints out a Java program that prints out the original program. Nothing new, an obvious generalisation of this. (Well, truth be told, it's a block of HTML that generates some Haskell code...)

But there's one big difference. To allow everyone else to join in the fun you can configure it yourself! Any non-empty list of languages will do, including length one. You may start hitting language line length limits if you make your list too long. Note that you can repeat languages so you can give someone some Haskell which decays into Perl after n iterations.

The code is geared towards generating tightly packed code but it's easily adapted to generate something slightly more readable. Apologies for the C warnings. Trivial to fix.

It's easily extended to support many more languages. C++, C#, obj-C, ocaml, Prolog, Lisp, Scheme and Go, say, should all be trivial apart from maybe a tiny bit of work with delimiting strings. (Eg. for Go you may need to tweak the import statement slightly so it doesnt't use double quotation marks.)

The code leaves many opportunities for refactoring but it's not like anyone is actually going to use this code for real production so I'm leaving it as is now.

I've only tested it under MacOS X. I don't know if there are carriage return/linefeed issues with other OSes. The shell script at the end is a regression test.

There was a little bit of theory involved which I learnt from Vicious Circles.

Here's a challenge for you: write a quine that takes as input the name of a language and outputs the same thing implemented in the input language. Much harder than what I just wrote.

> import Data.List

Here's the bit you can easily play with:

> langs = [Haskell, Perl, Python, Ruby, C, Java]

Implementation

> data Languages = Haskell | Ruby | Perl | C | Python | Java

> sequenceFromString Haskell s = "map toEnum[" ++ (intercalate "," $
>     map (\c -> show (fromEnum c)) s) ++ "]"
> sequenceFromString Perl s    = (intercalate "," $
>     map (\c -> "chr(" ++ show (fromEnum c) ++ ")") s)
> sequenceFromString Python s  = (intercalate "+" $
>     map (\c -> "chr(" ++ show (fromEnum c) ++ ")") s)
> sequenceFromString Ruby s    = (intercalate "+" $
>     map (\c -> show (fromEnum c) ++ ".chr") s)
> sequenceFromString C s       = concatMap
>     (\c -> "putchar(" ++ show (fromEnum c) ++ ");") s
> sequenceFromString Java s    = concatMap
>     (\c -> "o.write(" ++ show (fromEnum c) ++ ");") s

> paramList' Haskell = intercalate " " . map (\n -> "a" ++ show n)
> paramList' C       = intercalate "," . map (\n -> "char *a" ++ show n)
> paramList' Python  = intercalate "," . map (\n -> "a" ++ show n)
> paramList' Ruby    = intercalate "," . map (\n -> "a" ++ show n)
> paramList' Java    = intercalate "," . map (\n -> "String a" ++ show n)

> paramList Perl    _ = ""
> paramList lang n = paramList' lang [0..n-1]

> driver l args = defn l ++ intercalate (divider l) args ++ endDefn l

> divider C       = "\",\""
> divider Perl    = "','"
> divider Ruby    = "\",\""
> divider Python  = "\",\""
> divider Haskell = "\" \""
> divider Java    = "\",\""

> defn C       = "main(){q(\""
> defn Perl    = "&q('"
> defn Python  = "q(\""
> defn Ruby    = "q(\""
> defn Haskell = "main=q \""
> defn Java    = "public static void main(String[]args){q(\""

> endDefn C       = "\");}"
> endDefn Perl    = "')"
> endDefn Python  = "\")"
> endDefn Ruby    = "\")"
> endDefn Haskell = "\""
> endDefn Java    = "\");}}"

> arg Haskell n = "a" ++ show n
> arg Perl n    = "$_[" ++ show n ++ "]"
> arg C n       = "printf(a" ++ show n ++ ");"
> arg Python n  = "a" ++ show n
> arg Ruby n    = "a" ++ show n
> arg Java n    = "o.print(a" ++ show n ++ ");"

> argDivide Haskell l = "++" ++ sequenceFromString Haskell (divider l) ++ "++"
> argDivide Perl l    = ","    ++ sequenceFromString Perl (divider l) ++ ","
> argDivide C l       = sequenceFromString C (divider l)
> argDivide Python l  = "+" ++ sequenceFromString Python (divider l) ++ "+"
> argDivide Ruby l    = "+" ++ sequenceFromString Ruby (divider l) ++ "+"
> argDivide Java l    = sequenceFromString Java (divider l)

> argList lang1 lang2 n = intercalate (argDivide lang1 lang2) $
>     map (arg lang1) ([1..n-1] ++ [0])

> fromTo Haskell l n = "q " ++ paramList Haskell n ++ "=putStrLn$a0++" ++
>     sequenceFromString Haskell ("\n" ++ defn l) ++ "++" ++
>     argList Haskell l n ++ "++" ++ sequenceFromString Haskell (endDefn l)
> fromTo Perl    l n = "sub q {" ++ "print $_[0]," ++
>     sequenceFromString Perl ("\n" ++ defn l) ++ "," ++ argList Perl l n ++ "," ++
>     sequenceFromString Perl (endDefn l ++ "\n") ++ "}"
> fromTo Python  l n = "def q(" ++ paramList Python n ++
>     "): print a0+" ++ sequenceFromString Python ("\n" ++ defn l) ++
>     "+" ++ argList Python l n ++ "+" ++ sequenceFromString Python (endDefn l)
> fromTo Ruby    l n = "def q(" ++ paramList Ruby n ++
>     ") print a0+" ++ sequenceFromString Ruby ("\n" ++ defn l) ++
>     "+" ++ argList Ruby l n ++ "+" ++ sequenceFromString Ruby (endDefn l ++ "\n") ++ " end"
> fromTo C       l n = "q(" ++ paramList C n ++ "){" ++ "printf(a0);" ++
>     sequenceFromString C ("\n" ++ defn l) ++ argList C l n ++
>     sequenceFromString C (endDefn l ++ "\n") ++ "}"
> fromTo Java    l n = "public class quine{public static void q(" ++
>     paramList Java n ++ "){java.io.PrintStream o=System.out;o.print(a0);" ++
>     sequenceFromString Java ("\n" ++ defn l) ++ argList Java l n ++
>     sequenceFromString Java (endDefn l ++ "\n") ++ "}"

> main = do
>     let n = length langs
>     let langs' = cycle langs
>     putStrLn $ fromTo (head langs') (head (tail langs')) n
>     putStrLn $ driver (head langs') $ zipWith (\lang1 lang2 -> fromTo lang1 lang2 n)
>         (take n (tail langs')) (tail (tail langs'))

Regression Test
Assuming this article is stored in quineCentral.lhs.

runghc quineCentral.lhs>1.hs
cat 1.hs
echo "---------------------------------"
runghc 1.hs>2.pl
cat 2.pl
echo "---------------------------------"
perl 2.pl>3.py
cat 3.py
echo "---------------------------------"
python 3.py>4.ruby
cat 4.ruby
echo "---------------------------------"
ruby 4.ruby>5.c
cat 5.c
echo "---------------------------------"
gcc -o 5 5.c
./5>quine.java
cat quine.java
echo "---------------------------------"
javac quine.java
java quine>7.hs
cat 7.hs
diff 1.hs 7.hs

Generalising Gödel's Theorem with Multiple Worlds. Part II.

2010-12-24T08:27:00.000-08:00

Introduction
Last time we looked at a method for testing whether propositions of provability logic were valid by looking at the consequences of propositions within nested collections of worlds. This lends itself naturally to an algorithm that we can implement in Haskell. The diagrams I used last time are a variant on what are known as tableaux. But tableaux can be used in a number of ways and so we need code that is suitably generalised. In the code that follows I've ensured that the core algorithm has an interface that is suitable for carrying out the four tasks I'll demand of it. This means that the code leans more towards practicality than mathematical elegance.

And I apologise in advance: this post is mostly a bunch of implementation details. We'll get back to some mathematics next time. Nonetheless, by the end of this post you'll have working code to test the validty of propositions of provability logic.

Central to tableau algorithms is pattern matching. I'd like the Haskell pattern matcher do much of the work, and to increase its flexibility I'll need this extension:

> {-# LANGUAGE ViewPatterns #-}

We'll need these libraries too:

> import Control.Applicative
> import Control.Arrow
> import Control.Monad
> import Data.Function
> import List
> import Maybe
> import Text.Show

Logical Propositions
Now we'll need a bunch of logical operators. The first three are constructors for a proposition type and the rest are sugar to make expressions look nicer:

> infixr 1 :->
> infixr 2 :\/
> infixr 3 :/\

> infixr 1 -->
> infixr 1 <--
> infixr 1 <->
> infixr 2 \/
> infixr 3 /\

> (\/)    = (:\/)
> (/\)    = (:/\)
> (-->)   = (:->)
> (<--)   = flip (:->)
> p <-> q = (p :-> q) :/\ (q :-> p)

Here's our basic proposition type:

> data Prop = Letter String
>           | Prop :\/ Prop
>           | Prop :/\ Prop
>           | Prop :-> Prop
>           | Box Prop
>           | Dia Prop
>           | F
>           | T
>           | Neg Prop deriving (Eq, Ord)

I want show to know about operator precedence for propositions:

> instance Show Prop where
>     showsPrec p (a :/\ b) = showParen (p>3) $ showsPrec 3 a . showString " /\\ " . showsPrec 3 b
>     showsPrec p (a :\/ b) = showParen (p>2) $ showsPrec 2 a . showString " \\/ " . showsPrec 2 b
>     showsPrec p (a :-> b) = showParen (p>1) $ showsPrec 1 a . showString " --> " . showsPrec 1 b
>     showsPrec p (Neg r)   = showParen (p>4) $ showString "Neg " . showsPrec 5 r
>     showsPrec p (Box r)   = showParen (p>4) $ showString "Box " . showsPrec 5 r
>     showsPrec p (Dia r)   = showParen (p>4) $ showString "Dia " . showsPrec 5 r
>     showsPrec p (Letter n)= showParen (p>5) $ showsPrec 6 n
>     showsPrec p T         = showString "T"
>     showsPrec p F         = showString "F"

Some simple rules for simplification of some logical expressions:

> simplify p = let simplify' (a :\/ F) = a
>                  simplify' (F :\/ b) = b
>                  simplify' (a :/\ T) = a
>                  simplify' (T :/\ b) = b
>                  simplify' (a :\/ T) = T
>                  simplify' (T :\/ b) = T
>                  simplify' (a :/\ F) = F
>                  simplify' (F :/\ b) = F
>                  simplify' (F :-> b) = T
>                  simplify' (T :-> b) = b
>                  simplify' (a :-> F) = Neg a
>                  simplify' (a :-> T) = T
>                  simplify' (Neg T) = F
>                  simplify' (Neg F) = T
>                  simplify' (Box T) = T
>                  simplify' (Dia F) = F
>                  simplify' z = z
>    in case p of
>        a :/\ b -> let a' = simplify a
>                       b' = simplify b
>                   in simplify' (a' :/\ b')
>        a :\/ b -> let a' = simplify a
>                       b' = simplify b
>                   in simplify' (a' :\/ b')
>        a :-> b -> simplify' (simplify a :-> simplify b)
>        Box a   -> simplify' (Box (simplify a))
>        Dia a   -> simplify' (Dia (simplify a))
>        Neg (Neg a) -> simplify a
>        a           -> a

Kinds of Proposition
I'm actually going to use more than one proposition type. So when we do case analysis I need to make my patterns more abstract so they can work with multiple types. I'm going to use the PropType type to represent the ways we're going to classify logical propositions. The types are:

1. Atomic: A single letter or its negation
2. Constant: Simply T or F or a negation thereof.
3. DoubleNegation.
4. Disjunction: used to represent things like a∧b or ¬(a∨b).
5. Conjunction: used to represent things like a∨b or ¬(a∧b).
6. Provability: These are statements about provability like those starting with ◻ or ¬◊.
7. Consistency: These are statements about consistency like those starting with ◊ or ¬◻.

As this type is a simple container we can make it a functor too.

> data PropType a = Atomic a
>                 | Constant a
>                 | DoubleNegation a
>                 | Disjunction a a
>                 | Conjunction a a
>                 | Provability a
>                 | Consistency a

> instance Functor PropType where
>     fmap f (Atomic a)         = Atomic (f a)
>     fmap f (Constant a)       = Constant (f a)
>     fmap f (DoubleNegation a) = DoubleNegation (f a)
>     fmap f (Provability a)    = Provability (f a)
>     fmap f (Consistency a)    = Consistency (f a)
>     fmap f (Conjunction a b)  = Conjunction (f a) (f b)
>     fmap f (Disjunction a b)  = Disjunction (f a) (f b)

I'll introduce a typeclass that will allow us to use the PropType view to query propositions:

> class PropTypeable a where
>     propType :: a -> PropType a
>     neg      :: a -> a
>     isF      :: a -> Bool
>     negative :: a -> Bool
>     positiveComponent :: a -> Prop

And now here we have the cases that I summarised in English above:

> instance PropTypeable Prop where
>     propType (a :\/ b)        = Disjunction a b
>     propType (a :/\ b)        = Conjunction a b
>     propType (Neg (a :\/ b))  = Conjunction (Neg a) (Neg b)
>     propType (Neg (a :/\ b))  = Disjunction (Neg a) (Neg b)
>     propType (a :-> b)        = Disjunction (Neg a) b
>     propType (Neg (a :-> b))  = Conjunction a (Neg b)
>     propType (Neg (Neg a))    = DoubleNegation a
>     propType (Box a)          = Provability a
>     propType (Neg (Box a))    = Consistency (Neg a)
>     propType (Dia a)          = Consistency a
>     propType (Neg (Dia a))    = Provability (Neg a)
>     propType (Letter a)       = Atomic (Letter a)
>     propType (Neg (Letter a)) = Atomic (Neg (Letter a))
>     propType T                = Constant T
>     propType F                = Constant F
>     propType (Neg F)          = Constant T
>     propType (Neg T)          = Constant F
>     neg                       = Neg
>     isF F                     = True
>     isF (Neg T)               = True
>     isF _                     = False
>     positiveComponent (Neg a) = a
>     positiveComponent a       = a
>     negative (Neg _)          = True
>     negative _                = False

It'll be a while before we need the full generality so it's going to seem like overkill for the moment!

And some pre-packaged letters for convenience:

> [a, b, c, d, p, q, r, s, t] = map (Letter . return) "abcdpqrst"

We're going to need some operations that act on lists.

placesWhere finds all of the elements of a list for which some predicate holds. Instead of just listing the elements that match, it lists the elements paired with the rest of the list after the matching element is removed. We can think of these pairs as elements and their surrounding context:

> placesWhere p []     = []
> placesWhere p (x:xs) = let r = map (second (x:)) $ placesWhere p xs
>                        in if p x then ((x, xs) : r) else r

This finds something in the intersection of two sets using a given 'equality' predicate for matching. As we may be using a predicate different from == we need to see both of the (possibly different) elements that satisfy the predicate.

> findIntersection eq a b = listToMaybe [(x, y) | x <- a, y <- b, x `eq` y]

Sometimes we'll meet propositions that we can start to reason about using ordinary propositional calculus. These will match the propositional predicate:

> propositional (propType -> DoubleNegation _) = True
> propositional (propType -> Conjunction _ _)  = True
> propositional (propType -> Disjunction _ _)  = True
> propositional _                              = False

On the other hand the provability predicate is used to identify propositions that need rules pertaining to provability and consistency:

> provability (propType -> Provability _) = True
> provability (propType -> Consistency _) = True
> provability _                           = False

The Algorithm
And now we're almost ready to implement the tableau rules. Because we'll be using tableaux in a number of different ways I need lots of hooks into the algorithm that can perform different operations. I've collected all of these hooks into a single type. The algorithm will take a proposition of type prop (which will be Prop for the first three cases) and produce something of type result.

> data TableauRules prop result = TableauRules {

If our result corresponds to a world that is self-contradictory it is said to close. Here's how we indicate a closed (and hence not really existing) world. In the simplest case we won't actually store any information about a world, just whether or not it closes. So closes will be the identity function:

>     closes :: result -> Bool,

Occasionally we'll find a world with something obviously equivalent to F in it. It closes. Here's what we want to return in that case. The argument is the offending proposition:

>     foundF :: prop -> result,

Sometimes we'll find a pair that obviosuly contradict, like a and ¬a:

>     foundContradiction :: (prop, prop) -> result,

And sometimes we'll find an open world (ie. a real non-closed one). This function gets handed the list of propositions that have been found to hold in it:

>     open :: [prop] -> result,

Here's what we do when we find a conjunction. I hope you remember that when we meet a conjunction we can delete it and replace it with the two subpropositions. conjRule is handed the subpropositions as well as the result from proceeding with the tableau rule for conjunctions:

>     conjRule :: prop -> prop -> result -> result,

Disjunctions work a little differently. When handling a∨b, say, we need to handle two subtableaux, one with a and one with b. The first argument to disjRule rules is the disjunction itself, the next two are the left and right subpropositions, and the last two arguments are the results of continuing the two subtableaux.

>     disjRule :: prop -> prop -> prop -> result -> result -> result,

With doubleNegation we get to see propositions that have undergone double negation elimination.

>     doubleNegation :: prop -> result -> result,

When we use Dia to open new worlds we need to ensure that each of these subworlds is valid. Each subworld is processed with processWorld . For example, when we're drawing tableau diagrams we can use this book to draw a frame around the subtableaux. We then fold together these subworlds using combineWorlds:

>     processWorld  :: prop -> result -> result,
>     combineWorlds :: result -> result -> result,

Lastly we have our driver function that kicks off the whole tableau algorithm on a list of propositions:

>     tableau :: [prop] -> result -> result
> }

We can now implement a tableau computation algorithm that supports these hooks. We'll start with the check for whether or not we have an immediate contradiction:

> simpleClosure rules ps = case find isF ps of
>    Just a  -> foundF rules a

Split the propositions into those that are negated and those that are not. We're looking for propositions in one part that directly contradict propositions in the other:

>    Nothing ->
>      let (neg, pos) = partition negative ps
>          maybePair = findIntersection ((==) `on` positiveComponent) neg pos
>      in case maybePair of
>          Just pair -> foundContradiction rules pair
>          Nothing   -> open rules ps

Double negation elimination is straightforward to apply. Delete the original proposition and replace it with the version without double negation:

> applyDNeg rules p a props = doubleNegation rules a $
>   applyPropositional rules (a : delete p props)

The rule for handling conjunctions. We delete the conjunction from the current list of propositions and replace it with the two subpropositions:

> applyConj rules p a b props = conjRule rules a b $
>   applyPropositional rules (a : b : delete p props)

Disjunctions require running two separate subtableaux:

> applyDisj rules p a b props =
>    let props' = delete p props
>        left   = applyPropositional rules (a : props')
>        right  = applyPropositional rules (b : props')
>    in disjRule rules p a b left right

Here we tie together the rules for propositional calculus. We use a bit of case analysis to decide which rule to apply, and if no rule applies we try the provability rules instead:

> applyPropositional rules props =
>     let t = simpleClosure rules props in if closes rules t
>         then t
>         else case find propositional props of
>             Nothing -> applyProvability t rules props
>             Just p  -> case p of
>                 (propType -> DoubleNegation q) -> applyDNeg rules p q props
>                 (propType -> Conjunction a b)  -> applyConj rules p a b props
>                 (propType -> Disjunction a b)  -> applyDisj rules p a b props

When we've exhausted all possible rules from propositional calculus we scan for propositions like Dia p or Neg (Box p). These may imply the existence of subworlds. We then try to instantiate these subworlds, seeded according to the rules I gave in the previous article:

> applyProvability t rules props =
>     let impliedWorlds = placesWhere consistency props

>         consistency (propType -> Consistency _) = True
>         consistency _ = False

>         testWorld (p@(propType -> Consistency q), props) =

In the following line, neg p corresponds to the application of Löb's Theorem. provabilities is the list of propositions inherited by a subworld from statements about provability in the parent:

>              let tableau = runTableau rules (q : neg p : provabilities)
>                  provabilities = do 
>                      p@(propType -> Provability q) <- props
>                      [p, q]
>              in processWorld rules p tableau

>     in foldr (combineWorlds rules) t (map testWorld impliedWorlds)

And finally, here's where we kick our algorithm off:

> runTableau rules props = tableau rules props $ applyPropositional rules props

Testing Validity
We can now use a set of simple rules to test the validity of propositions. As mentioned above, the return data is a Bool used to indicate whether a subworld was invalid. By and large, these rules do the trivial thing. Note how the rule for disjunction requires both of the alternatives to be invalid in order to completely invalidate, so we use (&&). But when considering subworlds, just one bad subworlds is enough to invalidate a world:

> validRules = TableauRules {
>     closes = id,
>     open   = \_ -> False,

>     foundF             = \_ -> True,
>     foundContradiction = \_ -> True,

>     conjRule       = \_ _ t -> t,
>     disjRule       = \_ _ _ -> (&&),
>     doubleNegation = \_ t -> t,

>     combineWorlds = (||),
>     processWorld  = \_ t -> t,

>     tableau = \_ t -> t
> }

We can now write a simple validty test. We negate the proposition we're interested in and test whether the implied world closes:

> valid p = runTableau validRules [neg p]

Here's a small regression test to ensure everything works. It's just a bunch of examples that I worked out by hand or lifted from Boolos's book:

> valids = [
>         T,
>         a :-> a,
>         Box a :-> Box a,
>         Box a :-> Box (Box a),
>         Box (Box a :-> a) :-> Box a,
>         Box F <-> Box (Dia T),
>         let x = p :/\ q :-> r :-> a in Box (Box x :-> x) :-> Box x,
>         F :-> Dia p,
>         Box (Dia p) :-> Box (Box F :-> F),
>         (Box F \/ q /\ Dia (Box F /\ Neg q)) <->
>           (Dia (Box F \/ q /\ Dia (Box F /\ Neg q))
>           --> q /\ Neg (Box (Box F \/ q /\ Dia (Box F /\ Neg q)
>           --> q)))
>     ]

> invalids = [
>         F,
>         a :-> Box a,
>         Box a :-> a,
>         Box (Box a :-> a) :-> a,
>         Dia T,
>         Box (Dia T),
>         Neg (Box F),
>         (Box F \/ p /\ Dia (Box F /\ Neg q)) <->
>           (Dia (Box F \/ q /\ Dia (Box F /\ Neg q))
>           --> q /\ Neg (Box (Box F \/ q /\ Dia (Box F /\ Neg q)
>           --> q)))
>     ]

If everything is working, regress1 should give the result True.

> regress1 = do
>     print $ (and $ map valid valids) &&
>             (and $ map (not . valid) invalids)

That's enough implementation. In the next installment we'll start putting this code to work.

References
This code is an implementation of the algorithm described in Chapter 10 of The Logic of Provability. See that book if you want (1) a proof that the above algorithm always terminates and (2) that it really does correctly decide the validity of propositions of provability logic.

Generalising Gödel's Theorem with Multiple Worlds. Part I.

2010-12-11T15:58:00.000-08:00

Introduction
In his first incompleteness theorem, Gödel showed us that we can construct a sentence that denies its own provability. In his second incompleteness theorem he showed that an example of such a sentence is the one that asserts the inconsistency of arithmetic. If arithmetic is consistent then it can't prove its own consistency. On the other hand, if arithmetic is inconsistent then we can prove anything, and hence we can prove its consistency.

Can we generalise what Gödel did? For example, can we construct sentences that we can prove assert their own provability? What about sentences that deny that their provability is provable? Or what about sentences that assert that if they're provable then it's not provable that it's inconsistent that they imply that they're inconsistent with the rest of arithmetic?

Not only can we do these things, we can also write a computer program that generates such theorems for us. We can do so by working with the idea that a consistent set of axioms describes a world, and that a set of axioms able to talk about sets of axioms describes a world within a world, and a set of axioms that...well you can guess how it goes. A bit like Inception really.

Gödel's Theorem
Briefly, Gödel's first incompleteness theorem goes like this: we work with PA, the logical system built from Peano's Axioms. Within PA we can state and prove theorems about arithmetic like the fact that 1+2=3 or that for any prime there is always another greater prime. But even though PA is about numbers, we can talk about other things if we can encode them as numbers. In particular, Gödel came up with a scheme to encode propositions of PA as numbers. We use [P] to represent the number for the proposition P in Gödel's scheme. A proof is essentially just a list of propositions where each one is connected to some earlier ones by simple mechanical rules. These rules can be turned into arithmetical statements about the Gödel numbers for the propositions. So in the language of PA it is possible to assert that a particular number is the Gödel number of a provable proposition. Let Prov(n) denote the proposition of PA that says that n is the Gödel number of a provable proposition. Prov(n) is a proposition of PA and so is ¬Prov(n). Imagine we could find a proposition G with the property that G↔¬Prov([G]). It would assert its own unprovability. But it also appears to involve stuffing a representation of the Gödel number of G within G as well as all the rules for determining provability. Amazingly Gödel figured out how to do this (using tricks not unlike those used for quining). If G were false, G would be provable. So if we trust PA as a means of reasoning about proofs in PA, then G is true, though it can't be proved using the axioms of PA.

Provability Logic
We're going to be interested specifically in provability, so we don't need all of the power of PA. So we're going to work with a simplified domain specific logic called GL (for Gödel-Löb), otherwise known as Provability Logic. The idea is that sentences of GL will be shorthand for classes of statement in PA. GL will contain propositional calculus. So here's an example statement in GL: p∧q. The unknowns p and q represent propositions of PA and statements of GL are considered valid if they're provable whatever statements of PA they represent. For example, we could assign p="there are at least 10 primes" and q="7>2", in which case p∧q holds. But we could assign p="there are just 10 primes" in which case it's false. So p∧q isn't a valid proposition of GL as it doesn't hold for all p and q. On the other hand p→p∨q is valid because no matter what crazy propositions we assign to p and q, p→p∨q is true. (Of course, when I say "there are at least 10 primes", I mean a long and complicated sentence of PA that amounts to the same thing.)

But there's more to GL than propositional calculus. It also has the one-argument predicate ◻ which asserts that its argument is provable. More precisely, ◻p says that whatever proposition of PA p represents, let's say it's P, we have Prov([P]). It's just shorthand. Here's another example: ◻(q→◻p) says that whatever P we assign to p, and whatever Q we assign to q, Prov([Q →Prov([P])]). Or in English it says "for any propositions P and Q, it is provable that Q implies that P is provable".

We'll use ⊤ and ⊥ from propositional calculus as the always true and always false propositions in the usual way. ◻⊥ is the assertion that we can prove ⊥. In other words it's the assertion that PA is inconsistent. So now we can state Gödel's second incompleteness theorem: ¬◻⊥→¬◻¬◻⊥. If PA is consistent, we can't prove it.

We can also introduce the symbol ◊. This is just shorthand for ¬◻¬. ◊p says that it's not provable that p is false. In other words, it says that p is consistent with the rest of PA. So ◻ is provability and ◊ is consistency.

A set of assignments of propositions of PA to a bunch of letters in GL can be thought of as a world. For example, we can imagine a world in which p="2>1" and q="2<1". In this world, we have p and ¬q. The GL proposition ◻p says that we have a proof of p so it must be true in all worlds. Conversely, ◊p says that it's not true there is no world where p holds. In other words, it asserts the existence of a world where p holds. So worlds can talk about other worlds. If we have ◊◊p in some world, then it's asserting that there's another world in which ◊p holds. In other words, ◊◊p asserts there is a world in which it is asserted that there is another world where p holds. We can draw pictures to represent this. If a world has propositions that talk about another world then we draw the talked about world as a kind of subworld. Here's how we can picture ◊◊p:

Worlds are assignments of propositions of PA to letters of GL. But most of the time we won't particularly care about specific propositions themselves like "2>1". We'll be more interested in what truth assignments we can make to the propositions represented by the letters. So we can think of a world as a place where the letters of GL have been assigned truth values, true or false. And we can think of each world as containing subworlds consisting of propositions that can be proved or disproved by their parent worlds.

I'm going to spend most of my time looking at how we can explore and unfold all of the implications contained in a world.

The Rules of the Game
We can give some rules. If a world contains p∧q like this:

then it must contain both p and q as well:

I crossed out the p∧q as it became redundant. It's important to note that when I wrote p∧q we could have any propositions of GL standing in for p and q. So a world containing (p∨q)∧r must also contain p∨q and r. Other rules may apply too. If we had p∧q and r∧s we can unfold both to get p, q, r and s. This rule also kicks in if it applies to negated propositions that become conjunctions if we "push down" the negation using de Morgan's laws. For example if a world contains ¬(p→q) then it also contains p and ¬q.

If a world contains p∨q then we don't know exactly what it looks like.

There are two possibilities. We can draw both together like this:

The line in the middle means that the world either looks like what's on the left or it looks like what's on the right. Once we've indicated there are two alternatives then the original proposition became redundant again. When we have a vertical line like this then everything above the vertical line applies in both the left and right possibilities. This saves us having to split our world into two and copy all of our formulae to both sides. Like in the case of conjunctions this rule also kicks in for other kinds of disjunctions. So a world containing p→q splits into two separate worlds headed by ¬p and q.

If a world contains ⊥, or a contradiction, in means that it wasn't really a valid world after all. If we meet a world like this, we know that our starting point must have been been self-contradictory. If a world implies the existence of a subworld that isn't valid, then it must itself be invalid:

On the other hand, if we've split a world into two possibilities because we found p∨q then even if one branch is invalid the world might still be valid if the other branch is valid.

If we have a world containing ◻p:

then we know that p must be in every subworld. Actually, we know it must also hold in any subsubworld too, all the way down. This is because if we can prove something we can then use that proof directly to form a constructive proof that we can prove it. So that means ◻p holds in every subworld too:

We treat negation as above. So if we see ¬◻p we treat it like ◊¬p.

As we have also said, if we have a world with a ◊p in it then it must contain a subworld with p in it. Note that ◊ implies the existence of subworlds, but ◻ doesn't. It just tells us what must be in them if they exist.

But there's one last rule we'll need. It's Lob's theorem. This is the big theorem on which everything else I say depends. It states that ◻(◻p →p) →◻p. If we can prove that proving something implies its truth then we can prove it. I could sketch a proof here, but I highly recommend the cartoon proof here. In a way, Löb's theorem is a bit sad. The raison d'etre of mathematics is that we can use proofs to be sure of things. In other words, we take for granted that ◻p→p. But if we could prove this then we could prove ◻p, even if p were false! ◻p→p is not valid!

(Philosophical digression: Mathematicians assume ◻p→p. They don't assume ◻p→p because there's a proof. They assume it because experience has shown it to work. So I claim that ◻p→p is an empirical fact that we learn by scientific induction. That's controversial because I'm basically saying that much of mathematics is empirical - at least it is if you talk about truth rather than proof. If you disagree, don't worry about it. The rest of what I say here is independent of this digression.)

Anyway, we can flip Löb's theorem around to get ◊p →◊(p ∧¬◊p). So if we have a world in which ◊p holds:

then we have a subworld in which both p and ¬◊p hold:

And that's all we'll need apart from the obvious rule that we can remove double negation. Suppose we have a proposition of GL. We can use the rules above to extract as many implications as we can. Eventually there will come a point where there is nothing more we can do. If we can do this without hitting a contradiction then we've found a bunch of possible worlds for the proposition. If we can't find a valid world, however, the the original proposition must have been false. You might ask whether or not there any other rules we need in addition to the ones above. Maybe there are other theorems like Löb's theorem that we need to use. Amazingly it can be shown that no other rules are needed. This is a sure-fire terminating algorithm for determining whether or not a proposition of GL is valid! This is a powerful tool. We can now start constructing wild propositions like those I started with and find out whether they are valid. (Note that I've not in any way proved this procedure always works. You'll need to look at Boolos's book to see why.)

Some Proofs
Let's work through some examples. First I'll prove something from propositional calculus: (p∧q)∨(p∧r)→p. We start by drawing a world with the negation:

Now ¬(a→b) is the same as a∧¬b. So we get:

Now the only way to proceed is to consider the disjunction and consider two alternatives:

On both sides of the vertical line we find p. But that contradicts the ¬p we discovered earlier. So there is no way the negation of our original proposition can hold in any world. And therefore the original proposition must be valid.

How about showing Gödel's second incompleteness theorem this way. ◊⊤ says that ⊤ is consistent with the rest of arithmetic. Ie. it expresses the consistency of arithmetic. The theorem is then ◊⊤→¬◻◊⊤, ie. if PA is consistent, then PA can't prove it is consistent. We'll start with this:

We can deal with the negated implication like before. We can also deal straightforwardly with the double negation:

There's now no way to proceed except to use the ◊⊤ to open up a new world. Remember that when we use this we must use Löb's theorem as well as inheriting p and ◻p from the parent world's ◻p:

And we get a contradiction because we have both ◊⊤ and ¬◊⊤. So Gödel's second incompleteness theorem is indeed valid. Note that this isn't a demonstration from scratch. We've shown its validity from Löb's theorem. So this isn't really a useful way to show it's valid. But it *is* a useful way to show the validity of generalisations of Gödel's theorem. Unfortunately, I have to stop for now.

In the coming posts I'll implement all of the above as a computer program. If there's any ambiguity in what I've said I hope the source code to the program will resolve those ambiguities. What's more, our program won't just test the validity of a proposition but it will draw out a nice picture of our world with its subworlds. So you'll be able to trace through every step to make sure you understand!

But that's not all. We'll also see how to write a program to illustrate a bunch more theorems like Craig's Interpolation Lemma and Beth's Definability Theorem and then we'll finish with a program that is designed to construct self-referential propositions. In particular, given any self-referential description like "p is a proposition that is equivalent to the proposition that denies that it's provable that p's provability is consistent with p itself" it will solve to find p, even though GL doesn't allow us to directly construct self-referential propositions.

So let me recap: a world is an assignment of consistent truth values to the letters (and consequently propositions) of GL. Some of these propositions imply the existence of other subworlds with different truth values for these propositions. We draw these worlds as subworlds of the original world. For a world to be valid it mustn't contain any contradictions (unless the world contains a bunch of alternatives in which case just one alternative needs to be valid.) A proposition of GL is valid if unfolding its negation doesn't result in an invalid world.

Exercises
Which are the following are valid? If you report back any difficulties you have I can incorporate any needed revisions into the description above.
1. p→◻p
2. ◻p→◻◻p
3. ◊p→◊◻◻p
4. ◊(p→◻q) →◊(◻(◊(p∧q)))
5. ◻(p∧q∧r∧s) →◻p∧◻q∧◻r∧◻s

References
Most of this stuff is based on Chapter 10 of The Logic of Provability by Boolos. The idea of using multiple worlds to prove theorems is due to Kripke. I believe the procedure of unfolding the implications of a proposition in the tree-like way I describe above is due to Smullyan.

Update
Solutions to problems:
1. Not valid. Diagram:

2. Valid
3. Valid.
4. Valid.
5. Valid. Diagram:

Beating the odds with entangled qubits

2010-11-19T04:50:00.000-08:00

Quantum mechanics allows the possibility of "spooky action at a distance", correlations between widely separated but simultaneous random events that can't be explained by probability theory. These events look like they secretly communicate with each other, but we also know that quantum mechanics prevents us sending messages faster than the speed of light. Nonetheless, even though we can't exploit non-locality to send messages faster than the speed of light, two cooperating parties can exploit non-locality to perform tasks better than would be possible without non-locality. The CHSH game is one such example.

My goal here is to write code to emulate the CHSH game. It will require reusing the probability and quantum mechanics monads I've used here many times before. So I won't be explaining how these work.

> import Data.Map (toList, fromListWith)
> import Complex
> infixl 7 .*

The CHSH game is cooperative in the sense that the two players, A and B, are attempting to work together to win. The two players are widely separated. Between the two players is the game show host. The host randomly generates a pair of bits, s and t. s is sent to A and t is sent to B. Neither A nor B gets to see the message sent to their partner. A and B must now simultaneously make their moves, stating a choice of bit.

Call A's move a and B's move b. A and B win if they can arrange that a XOR b equals s AND t. So, for example, if A receives a true bit, and thinks B has also received a true bit, then A wants to make a move that differs from B's, otherwise A wants to make the same move as B. A and B are allowed to plan as much as they like beforehand but it should be pretty clear that they can't possibly guarantee a win.

We can formally write the victory condition as:

> victory a b s t = (a `xor` b) == (s && t)

Now there's a pretty good strategy A and B can adopt: three quarters of the time, s AND t will be false. In that case, A and B want their answers to match. So they could simply choose false, regardless of what message the game show host sends them. We can give their strategies as a function of the host's message:

> astrategy s = False
> bstrategy t = False

Now we can simulate our game:

> game = do

The host picks a random bit to send to each player:

>    s <- 0.5 .* return False + 0.5 .* return True
>    t <- 0.5 .* return False + 0.5 .* return True

The players now respond to each of their messages:

>    let a = astrategy s
>    let b = bstrategy t

Now we can collect the replies and score the result:

>    let score = victory a b s t
>    return score

> play1 = collect game

Running play1 gives the expected result that A and B have a 3/4 chance of winning. It's not hard to prove classically that they can do no better than this.

But in a quantum universe it is possible to do better! We now allow A and B to adopt strategies that involve making measurements of a quantum system. To describe their strategies we need to use the quantum monad. Here is the previous strategy rewritten for this monad. In this case, the argument b is the state of the quantum system they observe. The first element of the pair is the move in the game, the second element is the state the physical system is left in by the player:

> aqstrategy b s = qreturn (False, b)
> bqstrategy b t = qreturn (False, b)

Now we can rewrite game to support quantum processes:

> game' aqstrategy bqstrategy bits = do
>    s <- 0.5 .* preturn False + 0.5 .* preturn True
>    t <- 0.5 .* preturn False + 0.5 .* preturn True
>    (score, _, _) <- collect $ observe $ do
>        (abit, bbit) <- bits
>        (a, abit') <- aqstrategy abit s
>        (b, bbit') <- bqstrategy bbit t
>        let score = victory a b s t

Note that we have to return abit' and bbit' because quantum processes are reversible and can't erase information about a state. Also note that abit' and bbit' can be widely separated in space.

>        qreturn (score, abit', bbit')

Probabilistic processes can erase whatever they like:

>    preturn score

The quantum version of a coin toss to generate a random bit, ie. an equal superposition of False and True:

> coin' = (1/sqrt 2) .* qreturn False + (1/sqrt 2) .* qreturn True

> bits = do
>   a <- coin'
>   b <- coin'
>   return (a, b)

So now we play exactly as in play1 except that we give each player an independent qubit, which they ignore:

> play2 = collect $ game' aqstrategy bqstrategy bits

Unsurprisingly, the probability of winning is just the same as before.

But now we can try something impossible in the classical case. We give each player half of a perfectly correlated pair of qubits. The players are now entangled and can exploit non-locality. Of course if their strategies ignore the qubits we get the same result as before:

> bell = (1/sqrt 2) .* qreturn (False, False) + (1/sqrt 2) .* qreturn (True, True)
> play3 = collect $ game' aqstrategy bqstrategy bell

Now comes the surprising bit. The players can each look at the (classical) bit given to them by the game host. Depending what it is they rotate the qubit's state through some angle in state space. (In the case of the qubit being an electron state, this is an actual physical rotation of the electron.) In these strategies, the choice of move is the same as the state the qubit is left in after observation, hence the qreturn (b', b') bit.

> aqstrategy' b False = rotate (0) b >>= \b' -> qreturn (b', b')
> aqstrategy' b True  = rotate (pi/2) b >>= \b' -> qreturn (b', b')

> bqstrategy' b False = rotate (pi/4) b >>= \b' -> qreturn (b', b')
> bqstrategy' b True  = rotate (-pi/4) b >>= \b' -> qreturn (b', b')

And now when we play, the probability of winning is greater than 3/4.

> play4 = collect $ game' aqstrategy' bqstrategy' bell

All of the 'communication' took place before the game started. A and B didn't communicate s and t to each other. And yet they can beat the classical odds.

So in conclusion: Quantum mechanics gives opportunities for collusion that are impossible classically. Sadly we don't yet know how to maintain the state of separated entangled qubits for extended periods of time. But I remember seeing recently that people are managing to maintain qubit states for nanoseconds.

By the way, there are games where it is possible to achieve a 100% success rate with the help of quantum states. These give examples of what is known as quantum pseudo-telepathy. I presume the "pseudo" is because despite the 100% success rate, it still doesn't give a way to send messages instantly.

A last thought from me: one reason why humans send messages is to allow them to coordinate strategies. But quantum game theory shows that we can coordinate strategies without sending messages. In other words, even though non-locality doesn't give us faster-than-light communication, it does allow us to do things that were previously thought to require FTL. I think this may have some profound consequences.

And an example from a different domain: in biochemistry one could imagine remote parts of ligands coordinating the way they bind to receptors, something that would be completely missed by the kind of quasi-classical simulation I've seen biochemists use.

My standard quantum mechanics code:

> data W b a = W { runW :: [(a, b)] } deriving (Eq, Show, Ord)

> mapW f (W l) = W $ map (\(a, b) -> (a, f b)) l

> instance Functor (W b) where
>  fmap f (W a) = W $ map (\(a, p) -> (f a, p)) a

> instance Num b => Monad (W b) where
>  return x = W [(x, 1)]
>  l >>= f = W $ concatMap (\(W d, p) -> map (\(x, q)->(x, p*q)) d) (runW $ fmap f l)

> a .* b = mapW (a*) b

> instance (Eq a, Show a, Num b) => Num (W b a) where
>   W a + W b = W $ (a ++ b)
>   a - b = a + (-1) .* b
>   _ * _ = error "Num is annoying"
>   abs _ = error "Num is annoying"
>   signum _ = error "Num is annoying"
>   fromInteger a = if a==0 then W [] else error "fromInteger can only take zero argument"

> collect :: (Ord a, Num b) => W b a -> W b a
> collect = W . filter ((/= 0) . snd) . toList . fromListWith (+) . runW

> type P a = W Double a
> type Q a = W (Complex Double) a

> a `xor` b = a/=b

> rotate :: Double -> Bool -> Q Bool
> rotate theta True = let theta' = theta :+ 0
>   in cos (theta'/2) .* return True - sin (theta'/2) .* return False
> rotate theta False = let theta' = theta :+ 0
>   in cos (theta'/2) .* return False + sin (theta'/2) .* return True

> observe :: Ord a => Q a -> P a
> observe = W . map (\(a, w) -> (a, magnitude (w*w))) . runW . collect

Some help for the compiler (and maybe humans too):

> preturn = return :: a -> P a
> qreturn = return :: a -> Q a

Statistical Fingertrees

2010-11-09T07:26:00.000-08:00

I have no time to post a proper article. But I have time to post a mini-article with details in the links.

> {-# LANGUAGE MultiParamTypeClasses #-}
> import Data.FingerTree
> import Data.Monoid

I'm a bad person. I often use this method to compute the mean and variance of some samples. It's not robust. It performs badly if the variance is small compared to the size of the samples. I shouldn't use it. As penance, this is a quick article on dynamically computing variances robustly and efficiently for datasets that are frequently manipulated.

For convenience I'll talk about the unscaled variance: the sum, not the average of the square deviation from the mean. You can easily compute the variance from this if you know the dataset size.

If we know the size, mean and unscaled variance of two sets (more properly, multisets) we can find the mean and unscaled variance of the their union using the formula here. This method is much more robust that the naive algorithm.

The rule for combining two datasets gives us a monoid:

> data Stats = Stats { n :: Float, mean :: Float, unscaledVariance :: Float }
>              deriving Show

> instance Monoid Stats where
>     mempty = Stats 0 0 undefined
>     Stats n m v `mappend` Stats 0 _ _ = Stats n m v
>     Stats 0 _ _ `mappend` Stats n m v = Stats n m v
>     Stats n m v `mappend` Stats n' m' v' = 
>       let delta = m' - m
>       in Stats (n + n') ((n*m+n'*m')/(n+n'))
>                (v + v' + delta*delta*n*n'/(n+n'))

Given a single sample, we can compute its stats:

> instance Measured Stats Float where
>     measure x = Stats 1 x 0

Now we need just one more line of code:

> type StatsTree = FingerTree Stats Float

We now have a data structure that allows us to freely split, join, delete elements from and add elements to sequences of samples, all the while robustly keeping track of their mean and unscaled variance (and hence their variance).

For example:

> example = fromList [1..10] :: StatsTree
> test = let (_, b) = split ((>=4) . n) example
>            (c, _) = split ((>3)  . n) b
>        in measure c

computes the stats for the 3 elements starting at the 4th element of example.

An example application might be maintaining rolling averages and variances for a sliding window.

This was inspired an article by John D Cook somewhere around here.

On Removing Singularities from Rational Functions

2010-09-18T13:07:00.000-07:00

Introduction
Suppose we have the function

> f x = 1/(x+x^2) - 1/(x+2*x^2)

Some basic algebraic manipulation shows that in the limit as x→0, f(x)→1. But we can't simply compute f 0 because this computation involves division by zero at intermediate stages. How can we automate the process of computing the limit without implementing symbolic algebra?

I've already described one way to remove singularities from a function. But that approach is very limited in its applicability.

This article is about a variation on the approach to formal power series that nicely showcases some advantages of lazy lists. It will allow us to form Laurent series of functions so we can keep track of the singularities.

The usual Haskell approach to power series allows you to examine the coefficients of any term in the power series of the functions you can form. These series can't be used, however, to evaluate the function. Doing so requires summing an infinite series, but we can't do so reliably because no matter how many terms in a power series we add, we can never be sure that there aren't more large terms further downstream that we haven't reached yet. And if we want to perform computations completely over the rationals, say, we don't want to be dealing with infinite sums.

I'd like to look at a way of working with power series that allows us to perform exact computations making it possible to answer questions like "what is the sum of all the terms in this power series starting with the x^n term?" By extending to Laurent series, and implementing the ability to selectively sum over just the terms with non-negative powers, we can compute functions like f above at 0 and simply skip over the troublesome http://en.wikipedia.org/wiki/Pole_%28complex_analysis%29.

Power Series
When I previously discussed power series I used code that worked with the coefficients of the power series. This time we want to work with values of the function so it makes sense to store, not the coefficients a_i but the terms themselves, a_ixⁱ. So instead of a list of coefficients, Num a => [a] we need a representation that looks a little like:

> data Power a = Power (a -> [a])

where we pass x in as an argument to the function contained in a Power. But we also want to allow Laurent series so we need to also store an offset to say which (possibly negative) term our series starts with:

> data Laurent a = Laurent (a -> (Int, [a]))

But this fails us for at least two reasons:

1. We have the individual terms, but to evaluate the function requires summing all of the terms in an infinite list.
2. If we have a Laurent series, then we need to store values of a_ixⁱ for x=0 and i<0. We'll end up with division by zero errors.

Partial Sum Series
So here's what we'll do instead. Suppose our power series is Σ_i=n^∞a_ixⁱ. We'll store the terms s_j=Σ_i=j^∞a_ix^i-j. Our type will look like:

> data Partial a = Partial (a -> (Int, [a]))

> instance Eq (Partial a)
> instance Show (Partial a)

It's straightforward to add two functions in this form. We just add them term by term after first aligning them so that the xⁱ term in one is lined up with the xⁱ term in the other:

> instance Num a => Num (Partial a) where

>  Partial f + Partial g = Partial $ \x ->
>     let (m, xs) = f x
>         (n, ys) = g x
>         pad 0 _ ys = ys
>         pad n x ys = let z:zs = pad (n-1) x ys
>                      in x*z : z : zs
>         l = min m n
>     in (l, zipWith (+) (pad (m-l) x xs) (pad (n-l) x ys))

Notice the slight subtlety in the alignment routine pad. By the definition above, the jth term has a factor of x^j built into it. So we need to multiply by x each time we pad our list on the left.

Now we need to multiply series. We know from ordinary power series that we need some sort of convolution. But it looks like for this case we have an extra complication. We appear to need to difference our representation to get back the original terms, convolve, and then resum. Amazingly, we don't need to do this at all. We can convolve 'in place' so to speak.

Here's what an ordinary convolution looks like when we want to multiply the sequence of terms (a_i) by (b_i):

In this example, the blue diagonal corresponds to the terms that are summed to get the 4th term in the result.

However, we wish to work with partial sums s_j=Σ_i=j^∞a_ix^i-j and t_j=Σ_i=j^∞t_ix^i-j, constructing the partial sums of the convolution of a and b from s and t. The partial sums of the convolution can be derived from the partial sums by tweaking the convolution so it looks like this:

The blue terms work just like before and need to be summed. But we also need to subtract off the red terms, weighted by a factor of x. That's it! (I'll leave that as an exercise to prove. The inclusion-exclusion principle helps.)

The neat thing is that the red terms for each sum are a subset of the blue terms needed for the next element. We don't need to perform two separate sums. We can share much of the computation between the red and blue terms. All we need to do is write an ordinary convolution routine that additionally returns not just the blue terms, but a pair containing the blue sum and the red sum.

>  Partial f * Partial g = Partial $ \x ->
>     let (m, xs) = f x
>         (n, ys) = g x
>         (outer, inner) = convolve xs ys
>         f' a b = a-x*b -- (the subtraction I mentioned above)
>     in (m+n, zipWith f' outer inner)
>  fromInteger n = let n' = fromInteger n in Partial $ \_ -> (0, n' : repeat 0)

>  negate (Partial f) = Partial $ \x -> let (m, xs) = f x
>                                 in (m, map negate xs)
>  signum = error "signum not implemented"
>  abs    = error "signum not implemented"

This is an ordinary convolution routine tweaked to return the partial sum inner.

> convolve (a0:ar@(a1:as)) ~(b0:br@(b1:bs)) =
>  let (inner, _) = convolve ar br
>      ab = map (a0 *) bs
>      ba = map (* b0) as
>  in (a0*b0 : a0*b1+a1*b0
>            : zipWith3 (\a b c -> a+b+c) inner ab ba, 0 : inner)

The code is very similar to the usual power series multiplication routine. We can also use the same method described by McIlroy to divide our series.

As our series are a munged up version of the usual power series it's pretty surprising that it's possible to divide with so little code:

> instance Fractional a => Fractional (Partial a) where
>   fromRational n = let n' = fromRational n in Partial $ \_ -> (0, n' : repeat 0)
>   recip (Partial f) = Partial $ \x ->
>      let nibble (n, x:xs) | x==0      = nibble (n+1, xs)
>                           | otherwise = (n, (x:xs))
>          (n, xs) = nibble (f x)
>      in (-n, rconvolve x xs)

In effect, rconvolve solves the equation convolve a b==1:

> rconvolve x (a0:ar@(a1:as)) =
>   let (outer, inner) = convolve ar result
>       f a b = x*b-a
>       r = -1/f a0 a1
>       result = recip a0 : (map (r *) $ zipWith f outer inner)
>   in result

Note one ugly quirk of this code. I need to 'nibble' off leading zeroes from the series. This requires our underlying type a to have computable equality. (In principle we can work around this using parallel or).

That's it. We can now write a function to compute the positive part of a rational function. (By positive part, I mean all of the terms using non-negative powers of x.)

> pos f z = let Partial g = f $ Partial $ \x -> (1, 1 : repeat 0)
>               (n, xs) = g z
>           in if n>0
>               then z^n*head xs
>               else xs!!(-n)

Here are some examples:

> test1 = let f x = (1+2*x)/(3-4*x*x)
>         in pos (\x -> 1/(f x-f 0)/x) (0::Rational)

> test2 = pos (\x -> 1/(1+4*x+3*x^2+x^3) - 1/(1+x)) (1::Rational)

The original example I started with:

> test3 = pos (\x -> 1/(x+x^2) - 1/(x+2*x^2)) (0::Rational)

No division by zero anywhere!

Conclusions
The code works. But it does have limitations. As written it only supports rational functions. It's not hard to extend to square roots. (Try writing the code - it makes a nice exercise.) Unfortunately, any implementation of square root will (I think) require a division by x. This means that you'll be able to compute the positive part away from zero, but not at zero.

This method can't be extended fully to transcendental functions. But it is possible to add partial support for them. In fact, So with a little work we can still compute the positive part of functions like 1/sqrt(cos x-1) away from x==0. But applying cos to an arbitrary rational function may need more complex methods. I encourage you to experiment.

Note that this code makes good use of laziness. If your function has no singularities then you might find it performs no computations beyond what is required to compute the ordinary numerical value.

Automatic even/odd splitting

2010-09-05T09:34:00.000-07:00

Statement of the Problem
Suppose you have a real valued function on the reals, say f. We can split it into the sum of an even and odd part:

f(x) = f_odd(x)+f_even(x)

where

f_odd(x) = (f(x)-f(-x))/2, f_even(x) = (f(x)+f(-x))/2

If f_odd has a power series around zero, then all of its terms must have odd powers in x. So f_odd(x)/x must have all even powers and it becomes natural to 'compress' down the terms to form f'(x) = f_odd(sqrt(x))/sqrt(x). (Take that as a definition of f' for this article.) We can implement this operation as a higher order function:

> import Prelude hiding (odd)

> odd f x = let s = sqrt x in 0.5*(f s - f (-s))/s

Here's a simple example:

> f x = x*(1+x*(2+x*(3+7*x)))

> test0 = odd f 3

But there's something not quite right about this. If f is rational, then so is odd f. But the implementation of odd involves square roots. Among other things, square roots introduce inaccuracy. As square roots don't appear in the final result, can we eliminate them from the intermediate steps of the computation too?

A Better Solution
Let's use the results of last week's article to compute odd another way. We want a linear function that maps as follows:

1	->	0
x	->	1
x²	->	0
x³	->	x
x²ⁿ	->	0
x²ⁿ⁺¹	->	xⁿ

Here's an automaton:

If we start at 0, end at 1, and take exactly n steps, then the product of the factors we collect up along the way is given by the second column of that table. In n is even there is no such path so we collect up 0. As we're working with polynomials over the reals, rather than types, we have x1=1x and so on. We can construct a transition matrix:

0	1
x	0

We now do something like we did last time. Any time we have a function of some variable x, we replace x with the transition matrix. Our functions now take matrix values like

a	b
c	d

Any polynomial of our transition matrix always gives us equal elements along the diagonal. This is true even if we form the inverse of the transition matrix. So we don't need to store d. So now we implement a simple matrix type needing only to store three elements instead of four:

> data O a = O a a a deriving (Show, Eq)

> instance Num a => Num (O a) where
>    O a b c + O a' b' c' = O (a+a') (b+b') (c+c')
>    O a b c * O a' b' c' = O (a*a'+b*c') (a*b'+b*a') (c*a'+a*c')
>    fromInteger n = let i = fromInteger n in O i 0 0
>    negate (O a b c) = O (negate a) (negate b) (negate c)

Notice how similar this is to automatic differentiation. We can extend to reciprocals too:

> instance Fractional a => Fractional (O a) where
>    fromRational n = let i = fromRational n in O i 0 0
>    recip (O a b c) = let idet = recip (a*a-b*c) in O (idet*a) (idet*negate b) (idet*negate c)

And now we can implement a replacement for odd:

> transition x = O 0 1 x
> odd' f x = let O _ fx _ = f (transition x) in fx

> test1 = odd' f 3

Another example:

> test3 = odd' (\x -> (x+3*x*x-1/x)/(x*x)) 2
> test4 = odd  (\x -> (x+3*x*x-1/x)/(x*x)) 2

This new version has many advantages: it uses only rational functions, it's more accurate, and it's well defined at zero.

Conclusion
Automatic differentiation is just one of a family of methods that can be used to compute a wide variety of functions of real-valued functions. Essentially we're just working over real-valued matrices instead of real numbers. By using automata we can simplify the process of working out which matrices to use. (Though for the simple example above, you may have been able to guess the matrix without any other help).

(BTW I think there are hidden automata lurking in a few places in mathematics. For example, in Umbral calculus.)

Constraining Types with Regular Expressions

2010-08-14T17:31:00.000-07:00

Structures with Constraints
Here's a picture of an element of a binary tree type:

The leaves correspond to elements and the letters indicate the type of those elements. If we read the leaves from left to right they match the regular expression A^*B^*. We can define a binary type whose leaves always match this string by making use of dissections. Similarly we can construct all kinds of other structures, such as lists and more general trees, whose leaves match this expression.

But can we make structures that match other regular expressions? Here are some examples:

1. Structures that match the expression Aⁿ for some fixed n. These are structures, like trees, that have a fixed size. For the case n=5 we need a type that includes elements like:

2. Structures that match (AA)^*. These are structures like trees with an even number of elements. This is easy for lists. We can just use a list of pairs. But it's more complex for trees because subtrees may have odd or even size even though the overall structure is even in size. We could also generalise to ((AA)ⁿ)^* for some fixed n.

3. Structures that match A^*1A^*. Here's an example:

The element marked 1 is of unit type (). It's a 'hole'. So this regular expression is a derivative type. We can also construct a variety of types like those matching A^*BA^*CB^* by using a mixture of dissection and differentiation:

4. A bytestring- or rope-like type used to represent a string that is statically guaranteed to match a given regular expression.

5. Many kinds of other constraints you could imagine on a datastructure. Like trees which are guaranteed not to have two neighbouring leaves of the same type, or whose sequence of leaves never contain a certain subsequence.

So the challenge is this: can we implement a uniform way of building any container type to match any regular expression we want?

We'll need some extensions:

> {-# OPTIONS_GHC -fwarn-incomplete-patterns #-}
> {-# LANGUAGE TypeFamilies, EmptyDataDecls, UndecidableInstances, TypeOperators #-}

First Example: Differentiating Lists
Let l(x) be the type of lists of elements of type x. We can define a list type through

data List x = Nil | Cons x (List x)

Algebraically we can write this as

l(x) = 1+x l(x)

Let's try approaching the derivative of a list from the perspective of regular expressions. We know from Kleene's Theorem that the set of strings (ie. the language) matching a regular expression is precisely the language accepted by a finite state automaton. Let's consider the case x^*1x^*. This is the language we get by considering all paths from state 0 to state 1 in the following automaton:

Note that I have overloaded the numeral '1' to mean both the type with one element (when labelling an edge) and to mean the state numbered 1 (when labelling a vertex).

Define the L_ij(x) to be the type of lists whose sequence of elements correspond to all possible paths from state i to state j. A list in L_ij(x) is either an empty list, or constructed from a first element and a list. If we combine an element and a list, the combination has to match a possible path from i to j. There are a number of ways we could do this. The element could correspond to a transition from i to k. But if this is the case, then the remainder of the list must correspond to a path from k to j. So we must replace Cons with something that for its first argument takes a type corresponding to a single automaton step from i to k. For its second argument it must take an element of L_kj(x). The set of paths from i to j that correspond to a single step is the (i,j)-th element of the transition matrix for the automaton. This is X_ij(x) where

X(x) = (X_ij(x)) = (x 1)
(0 x)

On the other hand, we also need to replace Nil with a version that respects transitions too. As Nil takes no arguments, it must correspond to paths of length zero in the automaton. The only such paths are zero length paths from a state to itself. So the matrix for such paths is:

I = (1 0)
(0 1)

Let's also define the matrix L(x)= (L_ij(x)).

The words above boil down to:

L_ij(x) = I_ij+Σ_k X_ik(x)L_kj(x)

where the sum is over all the places we might visit on the first step of the journey from i to j.

We can rewrite this using standard matrix notation:

L(x) = I + X(x)L(x)

Compare with the definition of ordinary lists given above. We get the type of constrained lists by taking the original definition of a list and replacing everything with matrices. We replace 1 with I. We replace x with the transition matrix of the automaton. And we replace the structure we're trying to define with a family of types - one for each pair of start and end states for the automaton. We can describe this replacement more formally: it's a homomorphism from the set of types to the set of matrices of types. (Actually, it's a bit more subtle than that. This isn't quite the usual semiring of types. For one thing, the order of multiplication matters.) And it doesn't just apply to lists. We can apply this rule to any container type. For example, suppose we wish to repeat the above for trees. Then we know that for ordinary binary trees, t(x), we have

t(x) = x+t(x)²

We replace this with

T(x) = X(x)+T(x)²

We started this section by considering the specific case of the pattern x^*1x^* with a corresponding matrix X(x). Because X₁₀=0 and X₀₀=X₁₁ it's not hard to see that any type T we constrain using this regular expression will also have similar 'shape', ie. T₁₀=0 and T₀₀=T₁₁. So we can write

T = (t(x) t'(x))
(0 t(x))

where by definition, t(x)=T₀₀(x) and t'(x)=T₀₁(x). Suppose we have two such collections of types, (S_ij) and (T_ij). Now consider the types of pairs where the first element is of type S_ij and the second of T_jk. Then the leaves of the pair structure correspond to a path from i to k. So we have

(st(x) (st)'(x)) = (s(x) s'(x)) (t(x) t'(x))
(0 st(x)) (0 s(x)) (0 t(x))

Multiply out and we find that

(st)'(x) = s(x)t'(x)+s'(x)t(x)

In other words, the usual Leibniz rule for differentiation is nothing more than a statement about transitions for the automaton I drew above. To get a transition 0→1 you either go 0→0→1 or 0→1→1.

Although I talked specifically about differentiation, much of what I said above applies for any finite state automaton whose edges are labelled by types. The best thing now is probably to put together some code to see how this all looks.

A Specification
If you haven't checked Brent Yorgey's solution to my problem last week, now is a good opportunity. My code is a generalisation of that but it may be helpful to look at Brent's specialisation first.

The goal is to be able to define a transition matrix like this. ('K' is an abbreviation for 'constant' matrix.)

> type D a = K22 (a,    ())
>                (Void, a )

And then define a functor like this:

> type ListF = I :+: (X :*: Y)

Think of ListF being a bifunctor taking arguments X and Y.

We'd like then to able to form the matrix of fixed points of Y = ListF X Y. In this case, ordinary lists should appear as the element at position (0,0) in the matrix:

> type List x = Fix (I0, I0) (D x) ListF

I'm using In to represent the integer n at the type level.

Derivatives of lists should appear at (0,1) so we want

> type List' x = Fix (I0, I1) (D x) ListF

But Fix is intended to be completely generic. So it needs to be defined in a way that also works for trees:

> type TreeF = X :+: (Y :*: Y)
> type Tree  = Fix (I0, I0) (D Int) TreeF
> type Tree' = Fix (I0, I1) (D Int) TreeF

And of course it needs to work with other transition matrices. For example x^*1y^*1z^* has the following transition diagram and matrix:

> type E x y z = K33 (x,    (),   Void)
>                    (Void, y,    ()  )
>                    (Void, Void, z   )

So we'd expect

> type DDTree x y z = Fix (I0, I1) (E x y z) TreeF

to define the second divided difference of trees.

The Implementation

We'll need a type with no elements:

> data Void

And some type level integers:

> data Zero
> data S a
> type I0 = Zero
> type I1 = S I0
> type I2 = S I1
> type I3 = S I2

Now we'll need some type-level matrices. For any square matrix, we need a type-level function to give its dimension and another to access its (i, j)-th element:

> type family Dim x :: *
> type family (:!) x ij :: *

(Thanks to Brent's tutorial that code is much better than how it used to look.)

We can now define matrix addition through pointwise addition of the elements:

> data (:+) m n
> type instance Dim (m :+ n) = Dim m
> type instance (m :+ n) :! ij = Either (m :! ij) (n :! ij)

And similarly we can define multiplication. I'm using the type-level function Product' to perform the loop required in the definition of matrix multiplication:

> data (:*) m n
> type instance Dim (m :* n) = Dim m
> type instance (m :* n) :! ij = Product' I0 (Dim m) m n :! ij

> data Product' i k m n
> type instance Product' p I1 m n :! (i, j) = (m :! (i, p), n :! (p, j))
> type instance Product' p (S (S c)) m n :! (i, j) = Either
>    (m :! (i, p), n :! (p, j))
>    (Product' (S p) (S c) m n :! (i, j))
>

(Weird seeing all that familiar matrix multiplication code at the type level.)

Now we need some types to represent our functors:

> data I
> data X
> data Y
> data K n
> data (f :+: g)
> data (f :*: g)
> data F m f

(I think phantom empty data types should be called ethereal types.)

To turn these into usable types we need to implement the homomorphism I described above. So here are the rules laid out formally:

> type Id = K22 ((),   Void)
>               (Void, ())

> type family   Hom self m f :: *
> type instance Hom self m I = Id
> type instance Hom self m X = m
> type instance Hom self m Y = self
> type instance Hom self m (K n) = n
> type instance Hom self m (f :+: g) = Hom self m f :+ Hom self m g
> type instance Hom self m (f :*: g) = Hom self m f :* Hom self m g
> type instance Dim (F m f) = Dim m

> data Fix ij m f = Fix (Hom (F m f) m f :! ij)
> type instance (:!) (F m f) ij = Fix ij m f

That's more or less it. We can now go ahead and try to construct some elements. We could (as Brent suggests) write some smart constructors to make our life easier. But for now I'm writing everything explicitly so you can see what's going on:

> x0 = Fix (Left ())                         :: List Int
> x1 = Fix (Right (Left (1, Fix (Left ())))) :: List Int

x0 is the empty list. x1 is the list [1]. The Left and Right get a bit tedious to write. But this is intended as a proof that the concept works rather than a user-friendly API.

We can explicitly implement the isomorphism with the more familiar list type:

> iso1 :: [x] -> List x
> iso1 []     = Fix (Left ())
> iso1 (a:as) = Fix (Right (Left (a, iso1 as)))

> iso1' :: List x -> [x]
> iso1' (Fix (Left ()))              = []
> iso1' (Fix (Right (Left (a, as)))) = a : iso1' as
> iso1' (Fix (Right (Right a)))      = error "Can't be called as a is void"

So that's it! If we can write our container as the fixed point of a polynmomial functor, and if we can convert our regular expression to a finite state automaton, then Fix completely automatically builds the constrained container type.

What have we learnt?

We haven't just solved the original problem. We've shown that derivatives and dissections are special cases of a more general operation. Take a look at the definition of D x again. we can think of it as xI+Delta where Delta is the matrix

> type Delta = K22 (Void, ()  )
>                  (Void, Void)

This matrix has the property that its square is zero. It's the 'infinitesimal type' I described here. In other words, this is type-level automatic differentiation. We've also been doing type-level automatic divided differencing.

We can now go back and look at the matrix form of divided differences on wikipedia. I hope you can now see that the matrix T_id(x₀,...,x_n-1) defined there is nothing other than a transition matrix for this automaton:

In fact, we can use what we've learnt about regular expressions here to solve some numerical problems. But I won't write about that until the next article.

By the way, I think what I've described here can be viewed as an application of what Backhouse talks about in these slides.

I think that for any automaton we have a 2-category. The 0-cells are states, the 1-cells are the types associated with paths from one state to another, and the 2-cells are functions between types that respect the constraint. I haven't worked out the details however. The 2-category structure is probably important. As things stand, I've just shown how to make the types. But we don't yet have an easy way to write functions that respect these constraints. I suspect 2-categories give a language to talk about these things. But that's just speculation right now.

By the way, I couldn't write a working Show instance for Fix. Can you write one? And an implementation of arbitrary for QuickCheck?

And I hope you can now solve my problem from last week.

Leftover bits of code

K22 and K33 are constructors for 2×2 and 3×3 matrices. It would probably have been better to have used lists like Brent did.

> data K22 row0 row1
> type instance Dim (K22 row0 row1) = I2
> 
> type instance (:!) (K22 (m00, m01) row1) (I0, I0) = m00
> type instance (:!) (K22 (m00, m01) row1) (I0, I1) = m01
> type instance (:!) (K22 row0 (m10, m11)) (I1, I0) = m10
> type instance (:!) (K22 row0 (m10, m11)) (I1, I1) = m11
> 
> data K33 row0 row1 row2
> type instance Dim (K33 row0 row1 row2) = I3
> 
> type instance (:!) (K33 (m00, m01, m02) row1 row2) (I0, I0) = m00
> type instance (:!) (K33 (m00, m01, m02) row1 row2) (I0, I1) = m01
> type instance (:!) (K33 (m00, m01, m02) row1 row2) (I0, I2) = m02
> type instance (:!) (K33 row0 (m10, m11, m12) row2) (I1, I0) = m10
> type instance (:!) (K33 row0 (m10, m11, m12) row2) (I1, I1) = m11
> type instance (:!) (K33 row0 (m10, m11, m12) row2) (I1, I2) = m12
> type instance (:!) (K33 row0 row1 (m20, m21, m22)) (I2, I0) = m20
> type instance (:!) (K33 row0 row1 (m20, m21, m22)) (I2, I1) = m21
> type instance (:!) (K33 row0 row1 (m20, m21, m22)) (I2, I2) = m22

Update: I neglected to mention that there is a bit of subtlety with the issue of being able to create the same string by different walks through the automaton. I'll leave that as an exercise :-)

Divided Differences and the Tomography of Types

2010-08-07T15:52:00.000-07:00

Health Warning
This article assumes knowledge of Conor McBride's work on the differentiation and dissection of types, even if the early parts look deceptively simple.

Tables of Divided Differences
Given a sequence of numbers like 2, 6, 12, 20, 30 we can use the well known method of finite differences to predict the next in the series. We write the numbers in a column. In the next column we write the differences between the numbers in the previous column and iterate like so:

2
	4
6		2
	6		0
12		2
	8		0
20		2
	10
30

On the assumption that the rightmost column is zero all the way down we can extend this table to:

2
	4
6		2
	6		0
12		2
	8		0
20		2
	10		0
30		2
	12
42

We can think of this as an exercise in polynomial fitting. We have some x-values, in this case: x₀=1, x₂=2,..., x₄=5, and some y-values: y₀=2, y₁=6,..., y₄=30. We hope to fit a polynomial f so that f(x_i)=y_i. The table of finite differences has the property that the (n+1)-th column becomes constant if a degree n polynomial fits the data.

But what happens if we want to fit data to x-values that aren't equally spaced? Then we can use the *divided* difference table instead. In this case we don't record the differences between the y-values, but the quotient between the y-differences and the x-differences. For x_0,1,2,3 = 1, 2, 4, 6 and y_0,1,2,3 = 3, 6, 18, 138 we get:

1	3
		(6-3)/(2-1)=3
2	6		(6-3)/(4-1)=1
		(18-6)/(4-2)=6		(1-1)/(6-1)=0
4	18		(10-6)/(6-2)=1
		(38-18)/(6-4)=10
6	38

Again we reach zero. This is because y_i = f(x_i) = x_i²+2 so we have a quadratic again.

Let's assume y_i = f(x_i) for some set of points. We'll leave x, y and f unknown and fill out the table:

x₀	f(x₀)
		(f(x₁)-f(x₀))/(x₁-x₀)= f[x₀,x₁]
x₁	f(x₁)		(f[x₁,x₂]-f[x₀,x₁])/(x₂-x₀)= f[x₀,x₁,x₂]
		(f(x₂)-f(x₂))/(x₂-x₁)= f[x₁,x₂]		(f[x₁,x₂,x₃]-f[x₀,x₁,x₂])/(x₃-x₀)= f[x₀,x₁,x₂,x₃]
x₂	f(x₂)		(f[x₂,x₃]-f[x₁,x₂])/(x₃-x₁)= f[x₁,x₂,x₃]
		(f(x₃)-f(x₂))/(x₃-x₂)= f[x₂,x₃]
x₃	f(x₃)

Note that for any f, this table also defines the generalised divided differences f[x,y], f[x,y,z] and so on. (Compare with the notation at Wikipedia).

Divided Differences of Types
Now suppose that f isn't a function of real numbers, but is a container type. Then as I noted here, the second column corresponds to Conor McBride's dissection types. f[x,y] is an f-container with a hole, with everything to the left of the hole containing an x, and everything to the right of the hole containing a y.

So what's the next column? Consider f(x)=x⁴. Then f[x,y] = (y⁴-x⁴)/(y-x) = x³+x²y+xy²+y³. This corresponds to the description I just gave involving a hole. Now consider f[x,y,z] = (f[y,z]-f[x,y])/(z-x) = z²+xz+x²+y²+yz+yx. This is the type consisting of a 4-tuple, with two holes, and with x left of the first hole, y between the first and second holes, and z to the right of the second hole. In other words, it's a trisection of a container. More generally, f[x₀,...,x_n] is the type corresponding to n holes and blocks of elements of x_i between them. I've only shown this for the type x⁴ but it works for all the same types in Conor's paper. But there's a catch: I've used a definition involving subtraction, and that makes no sense for types. Don't worry, I'll address that later.

We can give trisections and so on a computational interpretation like in Conor's paper. Dissections correspond to tail recursion elimination. They make explicit the state of a function traversing and transforming a recursive type. (I recommend Conor's description in the paper so I won't reproduce it here.) Trisections correspond to a pair of coroutines. The first is transforming a recursive type. The second is transformed the result of the first routine. The second one can (optionally) start consuming pieces of tree as soon as the first has started producing them. At any moment in time we have 3 pieces: (1) the part of tree that is as yet untouched, (2) the part that the first coroutine has produced, but which the second hasn't seen yet, and (3) the output from the second coroutine. Make that explicit, and you get a trisection. The deal's much the same for quadrisections and so on.

Rediscovering Newton
A consequence of this is that we can now give a computational interpretation to much of the description of divided differences at Wikipedia. Among other things, when Conor derives a type theoretical analogue of the Taylor series, he's actually rediscovering a form of the Newton polynomial.

Type Matrices
One interesting property of divided differences is that they have a description in terms of matrices. In particular, Wikipedia describes a homomorphism from the ring of functions of x to the ring of matrices whose elements are functions of n unknowns. We know that types don't form a ring, but they do form a semiring. So we can form an analogous semiring homomorphism from the set of container types to the set of matrices whose entries are *types*. Matrices of types aren't something we see every day. How can we put this notion to work?

Let's use the notation from Wikipedia and consider the matrices we get when n = 1. (At this point I recommend reading all of the Wikipedia article on divided differences.) Let's define ourselves a tree container:

> data F a = Leaf a | Form (F a) (F a)

We'll call the homomorphism T so that

T(f)(x₀,x₁)=	f(x₀)	f[x₀,x₁]
	0	f(x₁)

Our tree type satisfies the equation f(x)=x+f(x)². As we have a homomorphism, we also expect this to hold:

T(f)(x₀,x₁) = T(i)(x₀,x₁)+T(f)(x₀,x₁)² (Equation 1)

where i is the identity function i(x) = x.

We can easily compute T(i)(x₀,x₁) directly using (y-x)/(y-x)=1. We get

T(i)(x₀,x₁)=	x₀	1
	0	x₁

Multiplying out the matrices in equation 1 now gives us a bunch of equations:

f(x₀) = x₀ + f(x₀)²
f(x₁) = x₁ + f(x₁)²
f[x₀,x₁] = 1 + f(x₀)f[x₀,x₁]+f[x₀,x₁]f(x₁)

The first two equations are just the definition of our tree. The third line, in Haskell, is:

> data F' x0 x1 = Empty | ForkL (F x0) (F' x0 x1) | ForkR (F' x0 x1) (F x1)

This is essentially the application of Conor's version of the Leibniz law to trees. F' is the dissected tree. And note how we don't have any subtraction. By using the matrix formulation of divided differences we only have matrix multiplications and sums to deal with.

The image of f(x) = x+f(x)² under the homomorphism T yields the simultaneous definitions of f and its dissection. More generally, if we had chosen to work with larger n we'd get the simultaneous definition of trees, their dissections, their trisections and so on. And it'll work for any recursive container type.

So matrices of types are a meaningful concept. They give a way to organise the mutually recursive definitions of higher order divided differences. If you look at this hard enough you may also see that the matrix I defined above is playing the same role as the D type in my previous article. With a little template Haskell I think we could in fact implement automatic (higher order) divided differences at the type level.

But this is all just scratching the surface. The matrix and homomorphism defined above don't just apply to divided differences. Matrices of types have another deeper and surprising interpretation that will allow me to unify just about everything I've ever said on automatic differentiation, divided differences, and derivatives of types as well as solve a wide class of problems relating to building data types with certain constraints on them. I'll leave that for my next article.

Exercise
In preparation for the next installment, here's a problem to think about: consider the tree type above. We can easily build trees whose elements are of type A or of type B. We just need f(A+B). We can scan this tree from left to right building a list of elements of type A+B, ie. whose types are each either A or B. How can we redefine the tree so that the compiler enforces the constraint that at no point in the list, the types of four elements in a row spell the word BABA? Start with a simpler problem, like enforcing the constraint that AA never appears.

Apology
Sorry about the layout above. My literate Haskell->HTML program doesn't support tables yet so there was a lot of manual HTML. This meant that I didn't write stuff out as fully as I could have and it may be a bit sketchy in places. I may have to switch to PDF for my next post. (Or use Wordpress, but I had trouble there too. or I could use one of the Javascript TeX renderers, but I don't like the external dependency.)

Automatic Divided Differences

2010-07-31T09:28:00.000-07:00

Divided Differences
I've previously talked about automatic differentiation here a few times. One of the standard arguments for using automatic differentiation is that it is more accurate than numeric differentiation implemented via divided differences. We can approximate f'(x) by using (f(x)-f(y))/(x-y) with a value of y near x. Accuracy requires y to be close to x, and that requires computing the difference between two numbers that are very close. But subtracting close numbers is itself a source of numerical error when working with finite precision. So you're doomed to error no matter how close you choose x and y to be.

However, the accuracy problem with computing divided differences can itself be fixed. In fact, we can adapt the methods behind automatic differentiation to work with divided differences too.

(This paragraph can be skipped. I just want to draw a parallel with what I said here. Firstly I need to correct the title of that article. I should have said it was about *divided differences*, not *finite differences*. The idea in that article was that the notion of a divided difference makes sense for types because for a large class of function you can define divided differences without using either differencing or division. You just need addition and multiplication. That's the same technique I'll be using here. I think it's neat to see the same trick being used in entirely different contexts.)

The Direct Approach
Firstly, here's a first attempt at divided differencing:

> diff0 f x y = (f x - f y)/(x - y)

We can try it on the function f:

> f x = (3*x+1/x)/(x-2/x)

diff0 f 1 1.000001 gives -14.0000350000029. Repeating the calculation with an arbitrary precision package (I used CReal) gives -14.000035000084000. We are getting nowhere near the precision we'd like when working with double precision floating point.

The Indirect Approach
Automatic differentiation used a bunch of properties of differentiation: linearity, the product rule and the chain rule. Similar rules hold for divided differences. First let me introduce some notation. If f is a function then I'll use f(x) for normal function application. But I'll use f[x,y] to mean the divided difference (f(x)-f(y))/(x-y). We have

(f+g)[x,y] = f[x,y]+g[x,y]
(fg)[x,y] = f(x)g[x,y]+f[x,y]g(y)
h[x,y] = f[g(x),g(y)]g[x,y] when h(x)=f(g(x))

We can modify the product rule to make it more symmetrical though it's not strictly necessary:

(fg)[x,y] = 0.5(f(x)+f(y))g[x,y]+0.5f[x,y] (g(x)+g(y))

(I got that from this paper by Kahan.)

In each case, given f evaluated at x and y, and its divided difference at [x, y], and the same for g, we can compute the corresponding quantities for the sum and product of f and g. So we can store f(x), f(y) and f[x,y] together in a single structure:

> data D a = D { fx :: a, fy :: a, fxy :: a } deriving (Eq, Show, Ord)

And now we can implement arithmetic on these structures using the rules above:

> instance Fractional a => Num (D a) where
>    fromInteger n = let m = fromInteger n in D m m 0
>    D fx fy fxy + D gx gy gxy = D (fx+gx) (fy+gy) (fxy+gxy)
>    D fx fy fxy * D gx gy gxy = D (fx*gx) (fy*gy) (0.5*(fxy*(gx+gy) + (fx+fy)*gxy))
>    negate (D fx fy fxy) = D (negate fx) (negate fy) (negate fxy)

I'll leave as an exercise the proof that this formula for division works:

> instance Fractional a => Fractional (D a) where
>    fromRational n = let m = fromRational n in D m m 0
>    D fx fy fxy / D gx gy gxy = D (fx/gx) (fy/gy) (0.5*(fxy*(gx+gy) - (fx+fy)*gxy)/(gx*gy))

For the identity function, i, we have i(x)=x, i(y)=y and i[x,y]=1. So for any x and y, the evaluation of the identity function at x, y and [x,y] is represented as D x y 1. To compute divided differences for any function f making use of addition, subtraction and division we need to simply apply f to D x y 1. We pick off the divided difference from the fxy element of the structure. Here's our replacement for diff0.

> diff1 f x y = fxy $ f (D x y 1)

This is all mimicking the construction for automatic differentiation.

Evaluating diff0 f 1 1.000001 gives -14.000035000083997. Much closer to the result derived using CReal. One neat thing about this is that we have a function that's well defined even in the limit as x tends to y. When we evaluate diff1 f 1 1 we get the derivative of f at 1.

I thought that this was a novel approach but I found it sketched at the end of this paper by Reps and Rall. (Though their sketch is a bit vague so it's not entirely clear what they intend.)

Both the Kahan paper and the Reps and Rall papers give some applications of computing divided diferences this way.

It's not clear how to deal with the standard transcendental functions. They have divided differences that are very complex compared to their derivatives.

Aside
There is a sense in which divided differences are uncomputable(!) and that what we've had to do is switch from an extensional description of functions to an intensional description to compute them. I'll write about this some day.

Note that the ideas here can be extended to higher order divided differences and that there are some really nice connections with type theory. I'll try to write about these too.

Update: I found another paper by Reps and Rall that uses precisely the method described here.

Death to Hydrae (or the operational semantics of ordinals)

2010-07-03T19:14:00.000-07:00

Unprovable Propositions
Among other things, Godel's first incompleteness theorem allows us to construct a statement in the language of Peano arithmetic that can't be proved using the axioms of Peano arithmetic. Unfortunately, this statement is a highly contrived proposition whose sole purpose is to be unprovable. People who learn of Godel's theorems often ask if there are other more natural and uncontrived mathematical statements that can't be proved from the Peano axioms.

My goal in this post will be to describe one of these propositions. Not just uncontrived, but actually very useful. I only intend to tell half of the story here because I feel like there are many good treatments already out there that tell the rest. I'm just going to get to the point where I can state the unprovable proposition, and then sketch how it can be proved if you allow yourself a little Set Theory.


> {-# OPTIONS_GHC -fno-warn-missing-methods #-}
> import Prelude hiding ((^))
> infixr 8 ^
> type Natural = Integer

Termination
Suppose we implement a function to compute the Fibonacci numbers like so:


> fib 0 = 0
> fib 1 = 1
> fib n = fib (n-2) + fib (n-1)

How do we know that fib terminates for all natural number arguments? One approach is this: if we pass in the argument n it clearly never recurses more than n levels. Each time it recurses it calls itself at most twice. So it must terminate in O(2ⁿ) steps (assuming that the primitive operations such as addition take constant time). We can think of this code in a kind of imperative way. It's a bit like n nested loops, each loop going round up to two times.

Suppose instead that we have some kind of recursive function g that goes n levels deep but for which the number of calls of g to itself is no longer two. In fact, suppose the number of self-calls is very large. Even worse, suppose that each time g is called, it calls itself many more times than it did previously, maybe keeping track of this ever growing number through a global variable. Or instead of a global variable, maybe an evil demon decides how many times g calls itself at each stage. Can you still be sure of termination?

A Simple Machine
In order to look at this question, we'll strip a computer right down to the bare minimum. It will have an input (that the evil demon could use) for natural numbers and will output only one symbol. Here's a design for such a machine:


> data Machine = Done | Output Machine | Input (Natural -> Machine)

A value of type Machine represents the state of the machine. Done means it has finished running. Output s means output a symbol and continue in state s. Input f means stop to input a number from the demon (or elsewhere), call it i, and then continue from state f i. This is very much in the style discussed by apfelmus and I in recent blog posts.

Here's an interpreter for one of these machines:


> run1 Done       = return ()
> run1 (Output x) = print "*" >> run1 x
> run1 (Input f)  = readLn >>= (run1 . f)

For any n we can easily build a machine to output n stars. This is such a natural machine to want to build it seems only right to give it the name n. If we want to do this then we need to make Machine an instance of Num and define fromInteger for it:


> instance Num Machine where
>    fromInteger 0 = Done
>    fromInteger n = Output (fromInteger (n-1))

Typing run1 8, say, will output 8 stars.

Now given two of these machines there is a natural notion of adding them. a + b is the machine that does everything b does followed by everything a does. (Remember, that's b then a.) To do this we need to dig into b and replace every occurrence of Done in it with a. That way, instead of finishing like b, it leads directly into a. In the case of a + Input f, for each number i we need to dig into f i replacing each Done with a:


>    a + Done     = a
>    a + Output b = Output (a + b)
>    a + Input f  = Input (\i -> a + f i)

There's a natural way to multiply these machines too. The idea is that in a * b we run machine b. But each time the Output command is run, instead of printing a star it executes a. You can think of this as a control structure. If n is a natural number then a * n means running machine a n times. In the case of a * Input f, instead of multiplying by a fixed natural number, we get an input from the user and multiply by f i instead:


>    _ * Done     = Done
>    a * Output b = a*b + a
>    a * Input f  = Input (\i -> a * f i)

We can make a machine to input a number and then output that many stars. Here it is:


> w = Input fromInteger

Try running run1 w.

Can you guess what the machine w * w does? Your first guess might be that it inputs two numbers and outputs as many stars as the product of the two numbers. Try it. What actually happens is that we're computing w * Input fromInteger. Immediately from the definition of * we get Input (\i -> w*i). In other words, the first input gives us an input i, and then w is run i times. So if we initially input i, we are then asked for i more inputs and after each input, the corresponding number of stars is output. Although the original expression contains just two occurrences of w, we are required to enter i+1 numbers.

Given the definitions of + and * it seems natural to define the power operation too:


> (^) :: Machine -> Machine -> Machine
> a ^ Done     = Output Done
> a ^ Output b = a^b * a
> a ^ Input f  = Input (\i -> a ^ f i)

The power operation corresponds to the nesting of loops. So, for example, w ^ n can be thought of loops nested n deep.

Try working out what w ^ w does when executed with run1.

Consider the set M of all machines built using just a finite number of applications of the three operators +, * and ^ to w and the non-zero naturals. (The non-zero condition means we exclude machines like 0*w that accept an input and do nothing with it.) Any such expression can be written as f w, where the definition of f makes no mention of w.

Suppose we use run1 and we always enter the same natural n. Then each occurrence of w acts like n. So if we start with some expression in w, say f w, then always inputting n results in f n stars. We could test this with run1 (w^w^w^w), always entering 2, but it would require a lot of typing. Instead we can write another intepreter that consumes its inputs from a list rather from the user (or demon). And instead of printing stars it simply prints out the total number of stars at the end:


> run2 Done    _        = 0
> run2 (Output x)   as  = 1 + run2 x as
> run2 (Input f) (a:as) = run2 (f a) as

Now you can try run2 (w^w^w^w) [2,2..] and see that we (eventually) get 2^{2^2²}.

Termination Again
If we run a machine in M there's a pattern that occurs again and again. We input a number, and then as a result we go into a loop requesting more numbers. These inputs may in turn request more inputs. Like the mythological hydra, every input we give may spawn many more requests for inputs. As the number of inputs required may depend on our previous inputs, and we may input numbers as large as we like, these machines may run for a long time. Suppose our machine terminates after requesting n inputs. Then there must be some highest number that we entered. Call it m. Then if the original machine was f w (with f defined in terms of the 3 operators and non-zero naturals), the machine must have terminated outputting no more than f m stars. So if our machine terminates, we can bound how many steps it took.

But do our machines always terminate? The input we give to the machine might not be bounded. If we run run2 (w^w) [4,5..], say, the inputs grow and grow. If these inputs grow faster than we can chop off the heads of our hydra, we might never reach termination.

Consider a program to input n and then output fib n. It accepts an input, recurses to a depth of at most n, and calls itself at most twice in each recursion. Compare with the machine 2 ^ w. This accepts an input n, recurses to a depth n, calling itself exactly twice each time. So if 2 ^ w terminates, so does fib. The more complex example above where I introduced the evil demon will terminate if w ^ w does, as long as the demon doesn't stop inputting numbers. So if we can show in one proof that every machine of type Machine terminates, then there are many programs whose termination we could easily prove.

Let's consider an example like run1 (w ^ w) with inputs 2, 3, 4, ...

We start with w ^ w. Examining the definition of the operator ^ we see that this proceeds by requesting an input. The first input is 2. Now we're left with w ^ 2. This is w * w. Again it accepts an input. This time 3. Now we go to state w * 3. This is w*2 + w. Again we accept an input. This time 4. We are led to w*2 + 4. This now outputs 4 stars and we are left with w * 2 which is w + w. We accept an input 5, output 5 stars and are left with w. After a further input of 6, it outputs 6 stars and terminates. Or we could just run run2 (w ^ w) [2,3..] and get 15(=4+5+6) as output.

The transitions are:
w ^ w
-> w ^ 2
-> w * w
-> w * 3
-> w*2 + w
-> w*2 + 4 -> ... -> w * 2
-> w + w
-> w + 5 -> ... w
-> 6 -> ... -> 0.

Now for some Set Theory. Rewrite the above sequence using the transfinite ordinal ω instead of w. The sequence becomes a sequence of ordinals. Any time we accept an input, the rightmost ω becomes a finite ordinal. So we have a descending sequence of ordinals. This is true whatever ordinal we start with. The execution of either Input or Output always strictly decreases our ordinal, and any descending sequence of ordinals must eventually terminate. Therefore every machine in M eventually terminates.

But here's the important fact: to show termination we used the ordinal ω, and this required the axiom of infinity and some Set Theory. Instead we could encode the termination question, via Godel numbering, as a proposition of Peano arithmetic. If we do this, then we hit against an amazing fact. It can't be proved using the axioms of Peano arithmetic. So we have here a useful fact, not a contrived self-referential one, that can't be proved with Peano arithmetic.

Why can't it be proved using just the Peano axioms?
A few years back, Jim Apple made a post about constructing (some) countable ordinals in Haskell. His construction nicely reflects the definitions a set theorist might make, but the code doesn't actually do anything. Later I learned from Hyland and Power how you can interpret algebraic structures as computational effects. apfelmus illustrates nicely how an abstract datatype can be made to do things with the help of an interpreter. Roughly speaking, doing this is what is known as operational semantics. So I thought, why not apply this approach to the algebraic rules for defining and combining ordinals. The result are the interpreters run1 and run2 above.

run1 gives an example of a Hydra game. In fact, its precisely the hydra game described in this paper because it always chops off the rightmost head. The Kirby-Paris theorem tells us we can't prove this game terminates using just the Peano axioms. A web search on Goodstein's theorem will reveal many great articles with the details.

A well-ordered quantity that you can keep decreasing as a program runs, and that can be used to prove termination, is an example of a loop variant. Loop variants are often natural numbers but the above shows that transfinite ordinals make fine loop variants. But in the interest of being fair and balanced, here's a dissenting view. The author has a point. If you are forced to use transfinite ordinals to show your program terminates, the age of the universe will probably be but the briefest flicker compared to your program's execution. On the other hand, if you don't want an actual bound on the execution time, ordinals can provide very short proofs of termination for useful programs.

(By the way, this article is a sequel to my own article.)

Exercise
Algebraic structures give rise to monads. Can you see how to generalise the definition of Machine to make it a monad? If you pick the right definition then the substitution of a for Done in the definition of + should give you a particularly simple definition of ordinal addition. (See this for a hint on how substitution works in monads.)


> instance Show Machine
> instance Eq Machine
> instance Ord Machine

Constructing Intermediate Values

2010-05-29T10:08:00.001-07:00

Introduction
The Intermediate Value Theorem, as it is usually told, tells us that if we have a continuous real-valued function f on the closed interval [0,1], with f(0)<0 and f(1)>0 then there is a point x in (0,1) such that f(x)=0. Here's a sketch of a proof:

Define three sequences l_i, u_i, x_i recursively as follows:

x_i+1 = (l_i+u_i)/2
(l_i+1, u_i+1) = (l_i, x_i+1) if f(x_i+1)>0
(l_i+1, u_i+1) = (x_i+1, u_i) otherwise

Set l₀=0, u₀=1.
The x_i+1 are bracketed by the l_i and u_i, and |u_i-l_i|=2^-i so the x_i form a Cauchy sequence with some limit x. Because f is continuous, we get f(x)=0.

This proof not only shows that we have an intermediate value x where f crosses the x-axis, it also seems to describe a procedure for computing x. Let's use this strategy to write a numerical method.

A Numerical Method

> import Ratio
> infixl 6 .*

> ivt0 :: (Float -> Float) -> Float
> ivt0 f = ivt' 0 1 where
>    ivt' l u =
>          let z = 0.5*(l+u)
>          in if z==l || z==u || f z == 0
>               then z
>               else if f z < 0
>                   then ivt' z u
>                   else ivt' l z

Here are some simple functions we can use to test it:

> f0, g0 :: Float -> Float
> f0 x = 3*x-1
> g0 x = 2*x-1

ivt0 f0 and ivt0 f1 give the expected results.

Of course, it might not actually be the case that f (ivt f) == 0. We're dealing with the type Float here and there might not actually be a solution to f x == 0 to the precision of a Float. And it's not clear that the word continuous means anything here anyway because the limited precision of a Float precludes us from constructing the deltas and epsilons the definition of continuity demands.

Switching to Exact Real Arithmetic
So can we do better? Can we implement a function to compute intermediate values for exact real arithmetic? I recently talked about computable real numbers, and they are what I intend to use here. I'll represent the real number x using a sequence of rationals, x_i, such that |x_i-x|<2^-i. Here's a suitable type:

> data R = R { runR :: Integer -> Rational }

(I'll only be calling it for non-negative arguments, but alas, Haskell has no natural number type.)

We'll need to display our numbers. Here's something cheap and cheerful that's good enough for this article:

> instance Show R where
>    show (R f) = show (fromRational (f 30) :: Float)

It displays a Float that is formed from a rational approximation that is within 2^-30 of our value.

We can inject the rationals into the computable reals by using constant sequences:

> instance Fractional R where
>    fromRational x = R $ const $ fromRational x

An important point here is that two different sequences might define the same real. So when we construct functions mapping reals to reals we must be sure that they take different representations of the same real to representations of the same real. I'll call these 'well-defined'.

Now our problems begin. Haskell wants us to implement Eq before Num, but equality of computable reals is not decidable. Still, it is semi-decidable in that if two reals differ, we can detect this. Here's a partial equality function:

> instance Eq R where
>    R f == R g = eq' 2 0 f g where
>        eq' delta i f g = 
>               let fi = f i
>                   gi = g i
>               in if fi<gi-delta
>                   then False
>                   else if fi>gi+delta
>                       then False
>                       else eq' (delta/2) (i+1) f g

The idea is that if f i and g i can be separated by more than 2*2^-i then the reals that f and g represent can't be equal.

We have to implement Ord too. If one computable real is bigger than another then they must differ by more than some power of 1/2. This means we can eventually distinguish between them. As a result, compare will always terminate if we compare distinct numbers:

> instance Ord R where
>    compare (R f) (R g) = compare' 2 0 f g where
>        compare' delta i f g =
>               let fi = f i
>                   gi = g i
>               in if fi<gi-delta
>                   then LT
>                   else if fi>gi+delta
>                       then GT
>                       else compare' (delta/2) (i+1) f g

We can define min and max straightforwardly:

>    min (R f) (R g) = R $ \i -> min (f i) (g i)
>    max (R f) (R g) = R $ \i -> max (f i) (g i)

And now we can define addition. We use the fact that if |x-x'|<ε and |y-y'|<ε then |(x+y)-(x'+y')|<2ε:

> instance Num R where
>    fromInteger n = R $ const $ fromInteger n
>    R f + R g = R $ \i -> let j = i+1 in f j + g j
>    negate (R f) = R $ negate . f

I'm going to skip multiplication of two reals. But it's easy to define multiplication of a rational with a real. Here I'm using the fact that if |x-x'|<ε, then |ax-ax'|<|a|ε:

> x .* R f = R $ mul x x f where
>    mul x x' f i | abs x' < 1 = x * f i
>                 | otherwise  = mul x (x'/2) f (i+1)

Now we're ready! It's straightforward to translate our previous algorithm. Note that we don't need a function to take the limit of a sequence. By definition we use sequences of rational approximations that are accurate to within 2^-i so our ivt1 function merely needs to compute an approximation to our intermediate value that is accurate to within 2^-i.

> ivt1 :: (R -> R) -> R
> ivt1 f = R $ \i -> ivt' 0 1 i where
>    ivt' :: Rational -> Rational -> Integer -> Rational
>    ivt' x y i =
>           let z = (1%2)*(x+y)
>           in if i==0
>                then z
>                else if f (fromRational z) < 0
>                    then ivt' z y (i-1)
>                    else ivt' x z (i-1)

Some new test functions:

> f, g :: R -> R
> f x = 3.*x-1
> g x = 2.*x-1

Note that (amazingly!) any well-defined function R -> R we implement is continuous because of the argument I sketched here.

And you'll find that ivt1 f gives the result you expect. But sadly, ivt1 g never returns. What went wrong?

Clutching Defeat from the Jaws of Victory
After one iteration, ivt1 sets z to the value we want. This seems like a particularly easy case. But the algorithm immediately fails because it then proceeds to compare f z with zero. We already know that our equality test must fail for this case. How do we fix our algorithm?

Consider the actions of any higher order function that operates on a continuous function f type R -> R. It will evaluate f at various real values. The return values from f will then be sampled at various values of i. Each sample corresponds to a bracket around a value of f. Consider a higher order function acting on the function below:

If our higher order function terminates, it can only compute a finite number of these brackets. So I can choose an ε small enough that I could nudge the graph up or down by ε without invalidating the brackets. I could provide the higher order function with a different function R -> R that responds in exactly the same way to the finite number of samples.

Now consider the case of an ivt function. It will respond the same way to two functions differing only by ε. The function in the graph above has been designed so that it is flat from 0.4 to 0.6. Nudging the function above by ε>0 up and down causes the crossing of the x-axis to move from x<0.4 to x>0.6. In other words, no algorithm could possibly return a valid intermediate value for all possible inputs because I can always contrive a pair of functions, to which ivt would respond identically, and yet which have widely separated points where they cross the x-axis. This shows that the problem of determining a crossing point of a continuous function that starts off less than zero, and ends up greater than zero, is uncomputable. We can do it for some functions. But no function can solve this problem for all continuous inputs. What's really weird about this is that it fails precisely when the problem seems easiest - when there are lots of zeroes!

A Revised Approach
But we don't need to give up. If we could ensure that at each bisection we could pick a point that wasn't exactly a zero, then we'd be fine. But how can we do this? Well here's a cheating way: we don't bother. We let that problem be the caller's responsibility. Instead of bisecting at each stage, we'll trisect. We'll pick a point that isn't a zero in the middle third by using a function pased in by the caller. Call this function g. g is of type Rational -> Rational -> Rational and the specification is that g x y must be a rational number between x and y that isn't a zero. Here's our new ivt2 function:

> ivt2 :: (Rational -> Rational -> Rational) -> (R -> R) -> R
> ivt2 g f = R $ \i ->
>     let delta = (1%2)^i
>         ivt' :: Rational -> Rational -> Rational
>         ivt' x y =
>               let l = (1%3)*(2*x+y)
>                   u = (1%3)*(x+2*y)
>                   z = g l u
>               in if u-l < delta
>                    then z
>                    else if f (fromRational z) < 0
>                        then ivt' z y
>                        else ivt' x z
>     in ivt' 0 1

How can we make use of this? Well suppose we want to find where a polynomial, f, crosses the x-axis. If the polynomial has degree n, it can only cross at n points. So if we pick n+1 values, f must evaluate to a non-zero value at one of them. We can examine the n+1 values of f(x) in parallel to find one that is non-zero. This must terminate. Here's an implementation for the case n+1=2. If we pass two real values to pickNonZero it will True or False according to whether the first or second argument is non-zero. It must only be called with a pair of numbers where one is non-zero:

> pickNonZero (R f) (R g) = pickNonZero' 0 1 where
>     pickNonZero' i e = let fi = f i
>                            gi = g i
>                        in if fi > e || fi < -e
>                           then True
>                           else if gi > e || gi < -e
>                               then False
>                               else pickNonZero' (i+1) ((1%2)*e)

Note how it examines the two values in parallel. If it simply tested the two numbers in turn against zero it would fail if the first were zero. (Interesting side note: we could write this using par.)

Suppose the function f is strictly monotonically increasing. Then we know that if y>x, then f(y)>f(x) so we can't have that both f(x) and f(y) are zero. We can use this to write a suitable helper function to write an intermediate value function that is 100% guaranteed to work for all monotonically increasing functions (assuming we don't run out of memory):

> monotonicPicker :: (R -> R) -> Rational -> Rational -> Rational
> monotonicPicker f x y = if pickNonZero
>                               (f (fromRational x)) (f (fromRational y))
>                             then x
>                             else y

> ivt3 :: (R -> R) -> R
> ivt3 f = ivt2 (monotonicPicker f) f

By varying monotonicPicker we can find crossing points for all manner of continuous function. Try ivt3 g. But no matter how hard we try, we can never write one that works for all continuous functions.

What went wrong?
So why did the Intermediate Value Theorem, despite its appearance of constructing an intermediate value, fail to yield an algorithm? For an answer to that, you'll have to wait until my next article. The surprising thing is that it's not a problem of analysis, but a problem of logic.

Acknowledgements
I have borrowed liberally from Paul Taylor, especially from here. However, I've chosen to phrase everything in terms of classical reasoning about Haskell programs. Any errors are, of course, mine. It was probably Andrej Bauer's awesome blog that got me interested in this stuff.

Afterword
I want to stress that I have made absolutely no mention of intuitionism or constructivism in what I've said above. I have simply used ordinary classical mathematics to reason about some Haskell programs. I'm leaving that for the next article.

Optimising pointer subtraction with 2-adic integers.

2010-05-15T11:49:00.001-07:00

Here is a simple C type and a function definition:


struct A
{
    char x[7];
};

int diff(struct A *a, struct A *b)
{
    return a-b;
}

It doesn't seem like there could be much to say about that. The A structure is 7 bytes long so the subtraction implicitly divides by 7. That's about it. But take a look at the assembly language generated when it's compiled with gcc:


movl 4(%esp), %eax
subl 8(%esp), %eax
imull $-1227133513, %eax, %eax
ret

Where is the division by 7? Instead we see multiplication by -1227133513. A good first guess is that maybe this strange constant is an approximate fixed point representation of 1/7. But this is a single multiplication with no shifting or bit field selection tricks. So how does this work? And what is -1227133513? Answering that question will lead us on a trip through some suprising and abstract mathematics. Among other things, we'll see how not only can you represent negative numbers as positive numbers in binary using twos complements, but that we can also represent fractions similarly in binary too.

But first, some history.

Some n-bit CPUs

That's an Intel 4004 microprocessor, the first microprocessor completely contained on a single chip. It was a 4 bit processor equipped with 4 bit registers. With 4 bits we can represent unsigned integers from 0 to 15 in a single register. But what happens if we want to represent larger integers?

Let's restrict ourselves to arithmetic operations using only addition, subtraction and multiplication and using one register per number. Then a curious thing happens. Take some numbers outside of the range 0 to 15 and store only the last 4 bits of each number in our registers. Now perform a sequence of additions, subtractions and multiplications. Obviously we usually expect to get the wrong result because if the final result is outside of our range we can't represent it in a single register. But the result we do get will have the last 4 bits of the correct result. This happens because in the three operations I listed, the value of a bit in a result doesn't depend on higher bits in the inputs. Information only propagates from low bit to high bit. We can think of a 4004 as allowing us to correctly resolve the last 4 bits of a result. From the perspective of a 4004, 1 and 17 and 33 all look like the same number. It doesn't have the power to distinguish them. But if we had a more powerful 8 bit processor like the 6502 in this machine, we could distinguish them:

Introduction
This is analogous to the situation we have with distances in physical world. With our eyes we can resolve details maybe down to 0.5mm. If we want to distinguish anything smaller we need more powerful equipment, like a magnifying class. When that fails we can get a microscope, or an electron microscope, or these days even an atomic force microscope. The more we pay, the smaller we can resolve. We can think of the cost of the equipment required to resolve two points as being a kind of measure of how close they are.

We have the same with computers. To resolve 1 and 17 we need an 8-bit machine. To resolve 1 and 65537 we need a 32-bit machine. And so on. So if we adopt a measure based on cost like in the previous paragraph, there is a sense in which 1 is close to 17, but 1 is even closer to 257, and it's closer still to 65537. We have this inverted notion of closeness where numbers separated by large (in the usual sense) powers of two are close in this new sense.

We have an interesting relationship between computing machines with different 'resolving' power. If we take an arithmetical computation on an N-bit machine, and then take the last M bits of the inputs and result, we get exactly what the M-bit machine would have computed. So an M-bit machine can be thought of as a kind of window onto the last M-bits onto an N-bit machine. Here's a sequence of machines:

Each machine provides a window onto the low bits of the previous machine in the sequence. But what happens at the "..." on the left? That suggests the bizarre idea that maybe all of these finite machines could be thought of as window to some infinite bit machine. Does that idea make any kind of sense?

I'll try to convince you that's a sensible idea by pointing out that it's something familiar to anyone who's taken a rigorous analysis course. (And I'll mention in passing that the above diagram illustrates a limit in an appropriate category!)

Mathematicians (often) build the real numbers from the rational numbers by a process known as completion. Consider a sequence like 1, 14/10, 141/100, 1414/1000, ... . The nth term is the largest fraction, with 10ⁿ in the denominator, such that its square is less than 2. It's well known that there is no rational number whose square is 2. And yet it feels like this sequence ought to be converging to something. It feels this way because successive terms in the sequence get as close to each other as you like. If you pick any ε there will be a term in the series, say x, with the property that later terms never deviate from x by more than ε. Such a sequence is called a Cauchy sequence. But these sequences don't all converge to rational numbers. A number like √2 is a gap. What are we to do?

Mathematicians fill the gap by defining a new type of number, the real number. These are by definition Cauchy sequences. Now every Cauchy sequence converges to a real number because, by definition, the real number it converges to is the sequence. For this to be anything more than sleight of hand we need to prove that we can do arithmetic with these sequences. But that's just a technical detail that can be found in any analysis book. So, for example, we can think of the sequence I gave above as actually being the square root of two. In fact, the decimal notation we use to write √2, 1.414213..., can be thought of as shorthand for the sequence (1, 14/10, 141/100, ...).

The notion of completeness depends on an an idea of closeness. I've described an alternative to the usual notion of closeness and so we can define an alternative notion of Cauchy sequence. We'll say that the sequence x₁, x₂, ... is a Cauchy sequence in the new sense if all the numbers from x_n onwards agree on their last n bits. (This isn't quite the usual definition but it'll do for here.) For example, 1, 3, 7, 15, 31, ... define a Cauchy sequence. We consider a Cauchy sequence equal to zero if x_n always has zero for its n lowest bits. So 2, 4, 8, 16, 32, ... is a representation of zero. We can add, subtract and multiply Cauchy sequences pointwise, so, for example, the product and sum of x_n and y_n has terms x_ny_n. Two Cauchy sequences are considered equal if their difference is zero. These numbers are called 2-adic integers.

Exercise: prove that if x is a 2-adic integer then x+0=x and that 0x=0.

There's another way of looking at 2-adic integers. They are infinite strings of binary digits, extending to the left. The last n digits are simply given by the last n digits of x_n. For example we can write 1, 3, 7, 31, ... as ...1111111. Amazingly we can add subtract and multiply these numbers using the obvious extensions of the usaul algorithms. Let's add ...1111111 to 1:


...11111111
...00000001
-----------
...00000000

We get a carry of 1 that ripples off to infinity and gives us zeroes all the way.

We can try doing long multiplication of ...111111 with itself. We get:


...11111111
...1111111 
...111111  
...11111   
...
-----------
...00000001

It's important to notice that even though there are an infinite number of rows and columns in that multiplication you only need to multiply and add a finite number of numbers to get any digit of the result. If you don't like that infinite arrangement you can instead compute the last n digits of the product by multiplying 11...n digits...111 by itself and taking the last n digits. The infinite long multiplication is really the same as doing this for all n and organising it in one big table.

So ...1111111 has many of the properties we expect of -1. Added to 1 we get zero and squaring it gives 1. It is -1 in the 2-adic integers. This gives us a new insight into twos complement arithmetic. The negative twos-complements are the truncated last n digits of the 2-adic representations of the negative integers. We should properly be thinking of twos-complement numbers as extending out to infinity on the left.

The field of analysis makes essential use of the notion of closeness with its δ and ε proofs. Many theorems from analysis carry over to the 2-adic integers. We find ourselves in a strange alternative number universe which is a sort of mix of analysis and number theory. In fact, people have even tried studying physics in p-adic universes. (p-adics are what you get when you repeat the above for base p numbers, but I don't want to talk about that now.) One consequence of analysis carrying over is that some of our intuitions about real numbers carry over to the 2-adics, even though some of our intuitive geometric pictures seem like they don't really apply. I'm going to concentrate on one example.

The Newton-Raphson Method
I hope everyone is familiar with the Newton-Raphson method for solving equations. If we wish to solve f(x)=0 we start with an estimate x_n. We find the tangent to y=f(x) at x=x_n. The tangent line is an approximation to the curve y=f(x) so we solve the easy problem of finding where the tangent line crosses the x-axis to get a new estimate x_n+1. This gives the formula

x_n+1 = x_n-f(x_n)/f'(x_n).

With luck the new estimate will be closer than the old one. We can do some analysis to get some sufficient conditions for convergence.

The surprise is this: the Newton-Raphson method often works very well for the 2-adic integers even though the geometric picture of lines crossing axes doesn't quite make sense. In fact, it often works much better than with real numbers allowing us to state very precise and easy to satisfy conditions for convergence.

Now let's consider the computation of reciprocals of real numbers. To find 1/a we wish to solve f(x)=0 where f(x)=1/x-a. Newton's method gives the iteration x_n+1 = x_n(2-ax_n). This is a well know iteration that is used internally by CPUs to compute reciprocals. But for it to work we need to start with a good estimate. The famous Pentium divide bug was a result of it using an incorrect lookup table to provide the first estimate. So let's say we want to find 1/7. We might start with an estimate like 0.1 and quickly get estimates 0.13, 0.1417, 0.142848, ... . It's converging to the familiar 0.142857...

But what happens if we start with a bad estimate like 1. We get the sequence:


1
-5
-185
-239945
-403015701065
-1136951587135200126341705

It's diverging badly. But now let's look at the binary:


00000000000000000000000000000000000000000000000000000000000000000000000000000000001
11111111111111111111111111111111111111111111111111111111111111111111111111111111011
11111111111111111111111111111111111111111111111111111111111111111111111111101000111
11111111111111111111111111111111111111111111111111111111111111111000101011010110111
11111111111111111111111111111111111111111111010001000101010011001000110110110110111
11100001111001111011011101100100000010000000011011010110110110110110110110110110111

Our series may be diverging rapidly in the usual sense, but amazingly it's converging rapidly in our new 2-adic sense!

If it's really converging to a meaningful reciprocal we'd expect that if we multiplied the last n digits of these numbers by 7 then we'd get something that agreed with the number 1 in the last 7 digits. Let's take the last 32 digits:


10110110110110110110110110110111

and multiply by 7:


10100000000000000000000000000000001

The last 32 bits are


00000000000000000000000000000001.

So if we're using 32 bit registers, and we're performing multiplication, addition and subtraction, then this number is, to all intents and purposes, a representation of 1/7. If we interpret as a twos complements number, then in decimal it's -1227133513. And that's the mysterious number gcc generated.

Epilogue
There are many things to follow up with. I'll try to be brief.

Try compiling C code with a struct of size 14. You'll notice some extra bit shifting going on. So far I've only defined the 2-adic integers. But to get he reciprocal of every non-zero number we need numbers whose digits don't just extend leftwards to infinity but also extend a finite number of steps to the right of the "binary point". These are the full 2-adic numbers as opposed to merely the 2-adic integers. That's how the extra shifts can be interpreted. Or more simply, if you need to divide by 14 you can divide by 2 first and then use the above method to divide by 7.

I don't know how gcc generates its approximate 2-adic reciprocals. Possibly it uses something based on the Euclidean GCD algorithm. I wasn't able to find the precise line of source in a reasonable time.

An example of a precise version of the Newton-Raphson method for the p-adics is the Hensel lemma.

The last thing I want to say is that all of the above is intended purely to whet your appetite and point out that a curious abstraction from number theory has an application to compiler writing. It's all non-rigourous and hand-wavey. Recommend reading further at Wikipedia. I learnt most of what I know on the subject from the first few chapters of Koblitz's book p-adic Numbers, p-adic Analysis, and Zeta functions. The proof of the von Staudt–Clausen theorem in that book is mindblowing. It reveals that the real numbers and the p-adic numbers are equally valid ways to approximately get a handle on rational numbers and that there are whole alternative p-adic universes out there inhabited by weird versions of familiar things like the Riemann zeta function.

(Oh, and please don't take the talk of CPUs too literally. I'm fully aware that you can represent big numbers even on a 4 bit CPU. But what what I say about a model of computation restricted to multiplication, addition and subtraction in single n-bit registers holds true.)

Some Challenges
1. Prove from first principles that the iteration for 1/7 converges. Can you prove how many digits it generates at a time?
2. Can you find a 32 bit square root of 7? Using the Newton-Raphson method? Any other number? Any problems?

Acknowledgements
Update: I have replaced the images with images that are, to the best of my knowledge, public domain, or composited from public domain images. (Thanks jisakujien. My bad.)

Update
If you want to play with some p-adics yourself, there is some code to be found here. That also has code for transcendental functions applied to p-adics.

Here's some C code to compute inverses of odd numbers modulo 2³² (assuming 32 bit ints). Like the real valued Newton method, it doubles the number of correct digits at each step so we only need 5 iterations to get 32 bits. (C code as I think it's traditional to twiddle one's bits in C rather than Haskell.)


#include <stdio.h>
#include <assert.h>

typedef unsigned int uint;

uint inverse(uint x)
{
   uint y = 2-x;
   y = y*(2-x*y);
   y = y*(2-x*y);
   y = y*(2-x*y);
   y = y*(2-x*y);

   return y;
}

int main()
{
   uint i;
   for (i = 1; i<0xfffffffe; i += 2)
   {
       assert (i*inverse(i) == 1);
   }
}

On representing some real numbers exactly

2010-04-24T08:47:00.000-07:00

There are uncountably many real numbers, and only countably many finite strings of symbols, so we know for sure that there is no scheme to represent all real numbers as finite strings of symbols in such a way that different reals get different representations. However, there are useful subsets of the real numbers that *are* equipped with finite encoding schemes. I hope to give some examples of such sets, including one, the set of periods, that isn't currently well known outside the mathematical world.

It's traditional to think of the real numbers as divided up into three types of number: the rational numbers, the algebraic numbers, and the leftovers known collectively as the transcendental numbers. If you're not familiar with those terms I'll be explaining them below. We can represent this as follows:

ℚ⊂A⊂ℝ

I'm using unicode for blackboard bold R and Q there. I hope your browser supports them. ℝ is the set of real numbers, ℚ is the set of rational numbers, and A is the algebraic numbers.

But there are many different kinds of transcendental number in existence, so the subject of this article will be three more levels of classification we could insert between the set of algebraic numbers and the set of real numbers in the sequence above.

But first, let's revisit ℚ and A.

The Rational Numbers, ℚ
The rational numbers are fairly straightforward to deal with. They are nothing more than the ratios of integers. Given two such ratios we can easily compare them to see whether they are equal. For example, we're taught at a young age how to tell if 2/3 is the same as 4/6.

However, it didn't take long for the Ancient Greeks to realise that they had some equations they could solve approximately using rationals, but couldn't solve exactly. The best known example is solving x²=2.

Assume x is rational, so that x=p/q, with p and q integers with highest factors cancelled out. We have p²=2q², and so p must be divisible by 2. But then p² is divisible by 4, and so q is divisible by 2. This means that we can cancel 2 from the top and bottom of p/q. That contradicts our assumption about cancellation. So the equation has no rational solution.

If we expect x²=2 to have a solution we must work with numbers that can't be represented as fractions.

Notice how the rational numbers come equipped with a notation to describe them. We can just write one integer over another, with a horizontal line between them. Given such a description we can immediately tell if it's valid (it's only invalid if the denominator is 0) and whether it equals some other rational number.

One interesting property of the rational numbers is that they are *dense* in the real numbers. This is just another way of saying that any real number can be approximated as well as we like by a rational number. For example, we can approximate √2 to within 1/1000 by the rational number 141421/10000. This is an important property that will come up again later.

We also know that there are only countably many rational numbers because there are only countably many pairs of integers.

The Real Algebraic Numbers, A
The real algebraic numbers are all the real numbers that we can obtain by finding roots of polynomial equations whose coefficients are rational numbers. For example, the real solutions to x⁵-2x²+½=0 are all algebraic numbers. Among other things, this includes the solution to x²=2 that wasn't contained in the rationals. The algebraic real numbers also include the solutions to "insoluble" algebraic equations. (By 'algebraic' equations, I just mean polynomial equations.) Even though Galois proved we can't write expressions for the solutions to all algebraic equations using the four arithmetic operations and nth roots, these solutions still exist in the real numbers, and they are all, by definition, algebraic numbers.

We can describe an algebraic number by writing down the equation it solves, and additionally providing some description to say which root of the equation we're interested in. Unfortunately, given any algebraic number there is always an infinite number of algebraic equations that it satisfies. So like with the rationals there is some redundancy. The good news, however, is that given two such descriptions of algebraic numbers there is an algorithm to tell whether or not they describe equal numbers.

Even though there are vastly more algebraic numbers than rational numbers, they still only form a countable set. As mentioned in the previous paragraph, every algebraic number can be described by a finite string of symbols, and there are only countably many such strings.

The Computable Real Numbers, C
Now I want to insert a set between A and ℝ. These are the computable real numbers. But first think back to what I said about the rational numbers.

Pick a real number x. Given any ε we can always find a rational number that comes within ε of x. We can define a function f on the rationals with the property that f(ε) is rational and within ε of x for any rational ε. So f gives rational approximations to x to any desired accuracy. Such a function uniquely specifies x. In fact x is the limit of f(ε) as ε tends to zero. The computable real numbers are the real numbers specified in this way by *computable* functions. Notice how there are no infinities involved here. A computable real number is represented by nothing more than a finite string of symbols forming a computer program that processes finite integers. We can find out what the real number is to any desired accuracy in a finite time.

Because we can write computer programs to find roots of algebraic equations to any desired accuracy we know that the algebraic numbers are contained in the computables. So now we have:

ℚ⊂A ⊂C ⊂ℝ

Where I'm using C to represent the computables.

Almost any number we ever want to work with is computable. Whether it's √2, or π, or the value of some humongous integral that's appeared in your engineering problem, chances are there's an algorithm somewhere that can approximate it as accurately as you like, given enough CPU time and RAM.

Note that every computable is described by a finite string of symbols, and so the computables form a countable set.

There are some problems with computable numbers. Given two representations of computable numbers we have no way of telling whether or not they are equal. It's not just that this is hard to do. There simply is no algorithm to prove that the numbers are equal. The problem is that the two computer programs may give exactly the same approximations to every degree of accuracy down to a certain ε. But for higher accuracy they may give different results. We have no way of knowing in advance at what value of ε they'll start differing. And anyway, there simply is no algorithm for telling if two computer programs will always generate the same results. So paradoxically, testing the equality of two computable numbers is itself uncomputable.

Actually, the situation is much worse. We can't even tell if we have a valid representation of a computable real. In order to do that, we need to know that our computer program to generate approximations terminates. But solving that would solve the halting problem.

By the way, despite these issues this representation of computable real numbers as functions leads to practical ways to do arbitrary precision arithmetic on a computer in a way that gives us guaranteed bounds on the accuracy of our results.

From here we could go two ways. We could try to rein in our class of number to something more manageable. Or we could try to push further and try to find ways to represent real numbers that aren't even computable. Let's do the former first.

The Periods, P
The periods are an interesting class of real number that doesn't seem to be well known. The definition may be old, but it was only recently that they were being promoted in mathematical circles as an interesting thing to study. So I thought I'd contribute to their promotion.

Consider the number π. It arises straightforwardly in geometric problems. An example is computing the area of a circle. Back in the 18th century Lambert proved that it was irrational and in the 19th century Lindemann proved that it was transcendental. But in a sense the transcendental numbers are simply a rubbish heap into which the leftover numbers have been discarded. Can we get some kind of handle on at least some of these numbers in a way that puts π back on the map?

The real number π is the area of a unit circle. So it is the area of the region of the plane given by the equation x²+y²<1.

Here's a similar construction of another real number:

The area of that region is log(2). (Natural logarithm of course!) We can construct the logarithm of any positive rational number in a similar way. These are all transcendental numbers. But in a way they are nice transcendental numbers. They arise from considering areas of the plane described by straightforward algebraic inequalities. (As in the previous section, I'm using 'algebraic inequalities' to mean 'polynomial inequalities'.)

That suggests a class of real number: those numbers that can be represented as the volume of a region defined by a bunch of algebraic inequalities with rational coefficients. These are known as the periods. I'm talking about generalised volumes in n-dimensions, not just areas in the plane. Clearly π and log(2) are periods.

But does this representation solve the problem we had with computable numbers where we were unable to guarantee we could check the equality of two numbers? Now it gets interesting. The answer is: we don't know. It is conjectured that if we have the same number represented in two different ways as a period that we can transform one representation into another simply by using a small set of elementary operations. It is also conjectured that there is a terminating algorithm for finding such a sequence of operations, or if one doesn't exist, demonstrating this.

One obvious question now is this: are there any computable reals that aren't periods? A recent paper, Periods and elementary real numbers claims to exhibit one by means of a kind of diagonalisation argument. But in general it's hard to prove that a number isn't a period. I don't believe there is a proof yet that e = 2.718... isn't a period, but mathematicians expect that it isn't.

I think the name 'period' comes from the fact that the period of a pendulum of rational length, in a gravitational field of rational strength, is a period. Computing this requires elliptic functions, and these often result in periods when applied to rational numbers.

Although mathematics is about numbers, few mathematical publications talk much about specific real numbers. Curiously, when they do, they often talk about numbers that are periods. For example there have been many recent papers on the values of the Riemann zeta function for interesting arguments. These are often periods.

It has been suggested that the study of periods is actually the study of algebraic geometry in disguise. There is certainly no end to the interesting mathematics we can do using only periods.

We now have:

Q ⊂A ⊂P ⊂C ⊂R

The Definables
In the section on computable numbers I mentioned that we could go the other way and try to find a bigger class of real numbers that could be represented by finite strings of symbols. We could simply go "all the way" and consider those real numbers that can be defined, by any means possible, using the symbols of mathematics. We can try to pin this down a bit better. We'll work with the language of set theory. In this language we can write strings of symbols like S = "x>0 and x²=2". Such a string uniquely defines a real number if when we glue the string "there exists a unique x such that" onto the beginning of it we get a true proposition. We can now represent the number √2 as the string S.

It looks like we now have:

Q ⊂A ⊂P ⊂C ⊂D ⊂R

with D the set of definables.

Except there's a problem. The set D is not definable! The problem is this: we can talk about strings of symbols representing real numbers. To talk about these in set theory we'd encode these strings of symbols as sets so that we can apply the language of set theory. In order to use our attempted definition of definable we need some way to say when a string of symbols represents a true proposition. But to define a set of definables we need to talk about true propositions in the language of set theory. Godel showed us how to talk about the provability of a proposition within set theory. He did this by showing that provability is about mechanical operations we can perform on strings. But there's nothing analogous for talking about the truth of propositions. In fact, Tarski showed us this is impossible. So while we can talk about all kinds of individual numbers as being definable, we can't construct the set of definable numbers.

You may be interested to see an example of a definable number that isn't computable. Probably the most publicised example is Chaitin's constant. It represents the probability that a randomly generated string of symbols is a computer program that terminates in a finite time. We can't actually compute this number because it requires us to solve the halting problem. Nonetheless, it's perfectly well defined.

You can find the original paper on periods by Zagier and Kontsevich here. I first found out about periods from the book Mathematics Unlimited.

Update: Jared asked if I intended the real algebraic numbers. That's what I said in the text, but I confusingly wrote ℚ for this set. But that usually means the *complex* algebraic numbers. So I've changed notation and now use A for the real algebraic numbers.

Similarly, the periods are usually defined to be complex numbers whose real and imaginary parts are given by algebraically specified volumes. I was trying to avoid mention of the complex numbers to keep prerequisites to a minimum and the essential ideas here work without mention of the complex numbers.

A Partial Ordering of some Category Theory applied to Haskell

2010-03-28T14:38:00.000-07:00

I've had a few requests from people wanting to teach themselves applications of Category Theory to Haskell based on posts I've made. I've made things difficult by posting stuff at random levels of difficulty and without any kind of organising thread through them. So here's an attempt to list a bunch of posts related to aspects of category theory. I've grouped them by themes and within each theme I've tried to list the articles in order of difficulty. Unfortunately there can be big gaps between one article and the next as none of this material was intended to be linked together continuously. Nonetheless, I hope this is of some help.

First a warning: Category Theory Screws You Up!.

Monads
The first theme has to be monads of course. But don't forget: monads don't do anything. They're simply an interface to something that you must already have implemented some other way. So don't believe all that hype about how monads are what allow Haskell to use side effects and I/O.

If you just want to use IO and don't care about monads: The IO Monad for People who Simply Don't Care

This gets more hits than anything else on my blog: You Could Have Invented Monads! (And Maybe You Already Have.)

The Trivial Monad

Homeland Security Threat Level Monad

Monads, Kleisli Arrows, Comonads and other Rambling Thoughts

The idea that monads give a way to describe substitutions in a tree forms the basis for this and the next two posts: Variable substitution gives a...

Monads are Trees with Grafting

Where do monads come from?

From Monoids to Monads

Some thoughts on reasoning and monads

The Mother of all Monads

Beyond Monads

Fold and Unfold
I think this is one of the great applications of Category Theory to Computer Science. Structural recursion can be characterised really nicely in terms of F-algebras. That's cool. But even cooler is that when you dualise the definitions you get a great way to look at non-terminating computations on things like streams.

Data and Codata

Coalgebras and Automata

Purely functional recursive types in Haskell and Python

Commutative Monads and Vector Spaces
Trying to order these is tricky. I'm not sure I define the term commutative monad until the talk.

I realised a while back that the operation to build a vector space from a basis is a monad. In fact, like many well known algebraic structures, we get a commutative monad. Monads for vector spaces, probability and quantum mechanics pt. I

Monads, Vector Spaces and Quantum Mechanics pt. II

Trace Diagrams with Monads

A talk this time: Commutative Monads, Diagrams and Knots

And what happens when you try to use a non-commutative monad when a commutative monad is expected: Why isn't ListT a monad?

Comonads
I'm not sure the killer application for comonads has been found yet. But I do think they're good for things like dataflow and cellular automata fit the comonad model very well:

Evaluating cellular automata is comonadic

Comonadic Arrays

Comonads and reading from the future

The Monads Hidden Behind Every Zipper

Categories of polynomials and comonadic plumbing

Category Theory
And these are generally categorical articles

You Could Have Defined Natural Transformations

Reverse Engineering Machines with the Yoneda Lemma

A Yonedic Addendum

Applying the Yoneda lemma to memoize functions that at first seem unmemoizable: Memoizing Polymorphic Functions with High School Algebra and Quantifiers

Products, Limits and Parametric Polymorphism

Dinatural Transformations and Coends

The goal here was not to understand what Category Theory applies to Haskell, but how Haskell code can be interpreted in other categories: "What Category do Haskell Types and Functions Live In?"

A first step towards 2-category theory: The Interchange Law

Products, Limits and Parametric Polymorphism

2010-03-06T09:26:00.000-08:00


> {-# LANGUAGE RankNTypes #-}

When I wrote about memoizing polymorphic types I mentioned that you can think of forall a. F(a) as the product over all types a of F(a), where F is some type level function. For example F might be a type constructor like []. That's not completely accurate, as I hope to now explain. Along the way we should get some insight into the meaning of the limit of a functor in a category.

Suppose we have two types, A and B. We can form their product (A, B). We have the two projections fst and snd and if we have an element x in (A, B) we know that there is no necessary relationship between fst x and snd x. We can freely choose x so that each of fst x and snd x can take on any values we like in A and B.

But now consider an element of forall a. F(a). For each concrete type X we have a projection π_X::forall a. F(a) -> F(X). So it looks like a product of all types. However, we can't freely choose elements of forall a. F(a) so as to get any element of F(X) we like for each choice of X. To demonstrate this, consider an element x of the type forall a. [a]. For any choice of X we get a projection. For example, picking X to be Int or String gives:


> p1 :: (forall a. [a]) -> [Int]
> p2 :: (forall a. [a]) -> [String]
> p1 x = x
> p2 x = x

We can draw a simple category diagram representing this:

However, forall a. [a] comes with a free theorem. For this particular type we have a free theorem that says that for any f :: X -> Y, fmap f (p1 x) == p2 x. For example consider the well known function show :: Int -> String. The free theorem tells us that this diagram commutes:

So if p1 x == [3] then p2 x == ["3"]. We have lost free choice. But we have lost a lot more freedom than this. We have a commuting triangle like this for absolutely any function f :: X -> Y. It should be clear that there is no way we can pick elements of our list to satisfy all of these constraints. So x must be the empty list.

Cones
This scenario of having one projection for each type has a name. It's an example of a cone. Let's borrow the definition from Wikipedia:

Let F:J→C be a functor. Let N be an object of C. A cone from N to F is a family of morphisms with one morphism for each X,

π_X:N→F(X)

so that for every morphism f:X→Y, the following diagram commutes:

For any Haskell functor F, the free theorem tells us that we have exactly these diagrams with forall a. F(a) playing the role of N. So forall a. F(a), with its projections to F(X), forms a cone.

Limits
If F is an instance of the Haskell Functor type class, ie. an endofunctor on Hask, then the type forall a. F(a) gives us a cone. But not just any old cone. I don't know how to prove this but I'm pretty sure it's true that the free theorems for a functor are the only non-trivial relations between F(X) and F(Y) that we're forced to obey. If that's true, then, in a sense, forall a. F(a) is the "biggest" type satisfying the free theorems. We can make this more precise by saying that for any cone N with associated projections π_X, we can can map it uniquely to forall a. F(a) so that the following diagram commutes:

This special kind of cone has a name. It's called a limit. In other words, ∀a. F(a) = lim F. In fact, this is exactly how Limit is defined in category-extras.

In a sense you can think of a limit of a functor in any category as being like a product for which a version of the free theorems for a functor holds.

Colimits
A dual story can be told for existential types and colimits. But to do this we need free theorems for existential types. I'll leave that until I've figured out a nice way to derive these free theorems...

Final words
I should have written this article before I wrote one on coends. Think of it as a prequel.

The fact that we don't have complete freedom of choice when defining polymorphic elements in Haskell is what we mean by 'parametric' polymorphism. Instead of specifying one individual value for each type we define the elements in a uniform way. In a language like C++ we can use template specialisation to freely construct a rule for getting a value from a type using 'ad hoc' polymorphism. That freedom comes at a price - it becomes harder to reason about polymorphism.

It's amazing that many definitions in category theory emerge naturally (pun fully intended) from the free theorems. I keep hoping that one day I'll find a paper on exactly what is going on here that I understand. The original free theorems paper is very uncategorical in its language.

(st(x)	(st)'(x))	=	(s(x)	s'(x))	(t(x)	t'(x))
(0	st(x))		(0	s(x))	(0	t(x))