How to write tolerably efficient optimzation code without really trying...

...or the (freely generated module over the (min,+)-semiring)-monad. (But if I used that as a title you probably wouldn't be reading this.)

There's a large class of optimisation problems that can be considered to be about graph searching. This can include anything from finding the shortest distance between two cities, through finding the edit distance between two strings to finding the most plausible reconstruction of a signal corrupted by noise. In each case there's an obvious implementation of a solution in Haskell: use the list monad to generate all of the paths in the search space and then find the optimal path. The problem with this is that generating all paths can result in combinatorial explosion with exponential running times. So this post is about what I think is a novel monad to fix that problem. It gives a simple DSL for specifying a wide range of optimisation problems with polynomial instead of exponential complexity. In some cases it results in something close to the best published algorithm with little effort on the part of the user.

This is quite a neat and useful monad. So if the following is hard going, skip to the examples first for some motivation. I'd hate people to miss out on something that turns out to be fairly easy to use once it's set up.

Unfortunately, despite the expressiveness of Haskell making this possible, we have to fight Haskell on two fronts. One is the old annoyance of Haskell not supporting restricted monads. For a solution to that I'll refer you to Eric Kidd who used an Oleg trick. So let's get the first one out of the way before talking about graphs. Here's a preamble that replaces the Prelude's monad class with our own:

> {-# LANGUAGE MultiParamTypeClasses,
>            NoImplicitPrelude,
>            UndecidableInstances #-}

> import Debug.Trace
> import Prelude hiding ((>>=),(>>),return)
> import Data.Monoid
> import Data.Map (toList,fromListWith)

> class Monad1 m a where
>  return :: a -> m a

> class (Monad1 m a, Monad1 m b) => Monad2 m a b where
>   (>>=) :: m a -> (a -> m b) -> m b
>   (>>) :: m a -> m b -> m b
>   l >> f = l >>= \_ -> f

> class MonadPlus m where
>    mzero :: m a
>    mplus :: m a -> m a -> m a

So here's a graph with lengths marked on the edges:

Suppose we want to find the shortest path from A to F in 3 hops. Before solving that, consider a more general problem. In the following graph

suppose we wish to find the shortest n-hop paths from A to the C_i. No matter what complexity is in the cloud, each of the paths must travel through one of the B_i. So we can solve the shortest paths problem by finding the shortest (n-1)-hop paths to the B_i and using:

dist_n(A,C_i) = min_j(dist_n-1(A,B_j)+d_ji)

Here's the crucial point to notice: If we consider dist_n(A,B_i) to be a vector whose indices are labelled by i, and d_ji to be a matrix, then the above formula is an ordinary matrix multiplication except that instead of addition and multiplication we use min and addition. If you've studied graph algorithms you'll recognise this as the same multiplication that appears in the Floyd-Warshall algorithm. If we extend the real numbers with infinity, then (R∪∞,min,+) forms a semiring. (See Cormen, Leiserson & Rivest for some good discussion of semirings in this context.) Just as you can form vector spaces over a field, you can form modules over this semiring. And just as vector spaces form a monad, so do these modules. But you don't need to know much about semirings because the monad itself is a fairly straightforward object.

Here's that infinity I was talking about:

> infinity :: Fractional a => a
> infinity = 1/0

I have to admit, I enjoyed writing that line of code :-) More importantly, Haskell does the right thing with infinities.

The monad is going to represent the length of a path to a graph vertex. Here's the definition:

> data M l p a = M { runM :: [(a,p,l)] } deriving (Eq,Show,Ord)

Ignore the l part for now, that's going to serve as a place to make a log, like the Writer monad. We'll need that when we want the actual shortest path, not just its length. For now we'll just use an empty list. So suppose that at some stage in our computation we've reached the node B with a shortest path of length 3 and node F with a shortest path of length 4. We represent this as

u = M [(B,3,[]),(F,4,[])] :: M Vertex Float [()]

.
Suppose we now consider what happens at the next step. We can express the options from B with

v = M [(A,1,[]),(D,3,[])]

and those from F as

w = M [(D,3,[]),(E,3,[])].

To combine these with the earlier computed paths we could replace B and F in u with v and w respectively to get:

M [(M [(A,1,[]),(D,3,[])],3,[]),(M [(D,3,[]),(E,3,[])],4,[])] :: M (M Vertex Float [()]) Float [()].

We need to reduce this back down to a M Vertex Float [()] using the matrix 'multiplication' formula above to pick the appropriate shortest path each time. But note the pattern: it's the usual M M a -> M a monad pattern. So we can write the substitution+reduction as an implementation of >>=. Here's the implementation, along with some other monadic goodness:

> instance (Ord a,Monoid l,Num p) => Monad1 (M l p) a where
>   return a = M [(a,0,mempty)]

> instance (Ord a,Monoid l,Fractional p,Ord b,Ord p,Ord l) => Monad2 (M l p) a b where
>   M x >>= f = let
>        applyf = [(f a,p,l) | (a,p,l) <- x]
>        rawresult = concat [[(a',(p+p',l `mappend` l')) | (a',p',l') <- a] | (M a,p,l) <- applyf]
>        reduced = map (\(a,(p,l)) -> (a,p,l)) $ toList $ fromListWith min $ rawresult
>        filtered = filter (\(_,p,l) -> p/=infinity) reduced
>        in M filtered

> instance MonadPlus (M l p) where
>    mzero = M mempty
>    mplus (M a) (M b) = M (a ++ b)

> guard True = return ()
> guard False = mzero

> choose x = M $ map (\(a,p) -> (a,p,mempty)) x

I've also had to provide my own implementation of MonadPlus. choose is a convenience function so we don't have to keep entering those empty lists.

So here's that graph represented in terms of a function that describes transitions from vertex to vertex:

> data Vertex = A | B | C | D | E | F deriving (Eq,Show,Ord)

> graph A = choose [(B,1),(C,2)]
> graph B = choose [(A,1),(D,3)]
> graph C = choose [(A,2),(D,4),(E,1)]
> graph D = choose [(B,3),(C,4),(E,2),(F,3)]
> graph E = choose [(C,1),(D,2),(F,3)]
> graph F = choose [(D,3),(E,3)]

And here's our first example. We compute all of the shortest 3-hop paths from A:

> ex0 :: M [()] Float Vertex
> ex0 = do
>    step1 <- return A
>    step2 <- graph step1
>    step3 <- graph step2
>    step4 <- graph step3
>    return step4

Just in case you think that's the wrong result note that these are the shortest paths to each accessible vertex that are exactly 3 hops long. If you want paths up to 3 hops long you have to add some zero length transitions from nodes to themselves. I'll leave that as a (trivial) exercise.

But suppose you just want the shortest path to F. Then you can use:

> ex1 :: M [()] Float Vertex
> ex1 = do
>    step1 <- return A
>    step2 <- graph step1
>    step3 <- graph step2
>    step4 <- graph step3
>    guard $ step4==F
>    return step4

Note how simple this code is. We just write code to walk the graph. All of the optimisation machinery is hidden in the background. But there's a dark secret here because the code isn't yet doing what a claim. But before I admit what the problem is, let's modify this code so we can actually see what the shortest paths are. record is just like tell from the Writer monad. I've implemented >>= so that it does the right thing:

> record l = M [((),0,l)]

> ex2 = do
>   step1 <- choose [(A,0),(B,0)]
>   record [step1]
>   step2 <- graph step1
>   record [step2]
>   step3 <- graph step2
>   record [step3]
>   step4 <- graph step3
>   record [step4]
>   guard (step4==E || step4==F)
>   return (step1,step4)

Just for fun I made this example more complex. It finds all four shortest paths from anywhere in {A,B} to anywhere in {E,F}. As an exercise try modifying the above code so that it finds the one shortest path from {A,B} to {E,F}. It's also simple to modify it to deal with conditions like the traveller being at one of a particular set of nodes at a particular time. It's a DSL, you have considerable freedom to form whatever conditions you want.

And now it's time to reveal the secret. Although I wrote >>= so that it reduces the list of considered paths at each stage the above code still has exponential complexity as you'll rapidly discover as you increase the number of steps.

Consider the third monad law: (m >>= f) >>= g = m >>= (\x -> f x >>= g).

If we evaluate the left hand side for this monad we get a reduction after applying f and then another after applying g. But we if we evaluate the right hand side we find that we get many small (and useless) reductions applied and then at the end a large reduction is applied. We don't get a big reduction at the intermediate step. The result: exponential complexity because the code ends up exploring every possible path performing stupid trivial reductions at each step. The monad law says we get the same result, but it doesn't take the same time.

And here's the really annoying thing: when you use do-notation in the usual way the Haskell compiler builds expressions like the slow right hand side. Bummer! But there's an easy, if tedious, fix. We have to use the third monad law as a rewrite rule to convert our code to the form on the left hand side. Here's ex2 rewritten:

> ex3 = do
>   step1 <- choose [(A,0),(B,0)]
>   step4 <- do
>       step3 <- do
>           step2 <- do
>               record [step1]
>               graph step1
>           record [step2]
>           graph step2
>       record [step3]
>       graph step3
>   do
>       record [step4]
>       guard (step4==E || step4==F)
>   return (step1,step4)

Not as pretty is it? It's still not too bad but it'd be nice to make this conversion automatic. I've a hunch this can be done using Haskell rules. Anyone care to volunteer some help making this work?

Anyway, I mentioned that lots of problems are in this class, some not obviously graph theoretical. Here's one: computing the edit distance between two strings:

> diag (a:as,b:bs) = choose [ ((as,bs),if a==b then 0 else 1) ]
> diag _           = mzero

> right (a:as,bs) = choose [((as,bs),1)]
> right _         = mzero

> down (as,b:bs) = choose [((as,bs),1)]
> down _         = mzero

> nop ([],[]) = return ([],[])
> nop _       = mzero

> munch :: (String,String) -> M [()] Float (String,String)
> munch x = diag x `mplus` right x `mplus` down x `mplus` nop x

> pure (M a) = length a==1
> editDistance a b = until pure (>>= munch) (munch (a,b))

> ex4 = editDistance "hypocrisy" "hippocracy"

I think this is actually the dynamic programming algorithm in disguise (and slightly reordered) with some extra log time penalties because it's using a binary tree (from Data.Set) instead of arrays. But I wouldn't stake my life on it.

Anyway, here's one last example. Rather than explain it I refer you to the nicely written Wikipedia on the Viterbi algorithm. In fact, I decided to lift their example precisely. And note that I've (1) written the code twice, the second one is the one rewritten for speed, (2) used logarithms to convert the multiplications to additions and (3) 'unrolled' the loop to make absolutely explicit what's going on. If I wrote Haskell code to work on general arrays like the Python code in the Wikipedia article, it'd be too short and you might miss it if you blinked.

> data Weather = Rainy | Sunny deriving (Eq,Show,Ord)
> data Activity = Walk | Shop | Clean deriving (Eq,Show,Ord)

> transition Rainy = choose [(Rainy,-log 0.7),(Sunny,-log 0.3)]
> transition Sunny = choose [(Rainy,-log 0.4),(Sunny,-log 0.6)]

> emission Rainy = choose [(Walk,-log 0.1),(Shop,-log 0.4),(Clean,-log 0.5)]
> emission Sunny = choose [(Walk,-log 0.6),(Shop,-log 0.3),(Clean,-log 0.1)]

> ex5 = do
>       day1 <- choose [(Rainy,-log 0.6),(Sunny,-log 0.4)]
>       record [day1]
>       emission1 <- emission day1
>       guard (emission1 == Walk)
>       day2 <- transition day1
>       record [day2]
>       emission2 <- emission day2
>       guard (emission2 == Shop)
>       day3 <- transition day2
>       record [day3]
>       emission3 <- emission day3
>       guard (emission3 == Clean)
>       day4 <- transition day3
>       record [day4]

> ex6 = do
>       day4 <- do
>           day3 <- do
>               day2 <- do
>                   day1 <- choose [(Rainy,-log 0.6),(Sunny,-log 0.4)]
>                   do
>                       record [day1]
>                       do
>                           emission1 <- emission day1
>                           guard (emission1 == Walk)
>                   transition day1
>               do
>                   record [day2]
>                   do
>                       emission2 <- emission day2
>                       guard (emission2 == Shop)
>               transition day2
>           do
>               record [day3]
>               do
>                   emission3 <- emission day3
>                   guard (emission3 == Clean)
>           transition day3
>       record [day4]

[Sunny,Rainy,Rainy,Rainy]. It's the Viterbi algorithm for free.

Anyway, I haven't proved a thing above. I think (hope?) I'm right about the polynomial complexity. Certainly on paper it all works, but I've always found it hard to reason about the complexity of Haskell code.

Unfortunately this code is all still slower than I'd like. I've played with trace during the development and everything seemed to run the way I expect (at least after I discovered the third monad law problem). I think that ultimately there's a big overhead for piling up HOFs, and the lookups in Data.Set probably aren't helping. Still, I don't think that invalidates the mathematics behind the approach and I think that asymptotically the slowness is really a big constant and a bunch of logs and it's still much better than combinatorial search. I ought to see how ocaml fares.

BTW There's a deep connection between this monad and the quantum mechanics one. I've no time to write about that now, though see here.

Some acknowledgements: the strange but highly readable Wikipedia article motivated me to think about the Viterbi algorithm, and I became sure there was a monad lurking in the Python code. I spent a month or so getting nowhere. And then I spoke with Eric Kidd on the phone about probability monads. He challenged me on something and in an attempt to formulate a response the whole of the above seemed to just crystallise. Thanks for asking the right question Eric!

I think there might be a comonad for working backwards through the graph...

Anyway, one last thing: please for the love of the FSM will someone put a clean way to implement restricted monads into their Haskell compiler. There is so much neat stuff that could be implemented with it and I now officially consider Haskell broken without one.

Haskell Incarnate

I'm off to the Make Faire in a few minutes and I'm running in Bay to Breakers tomorrow so I probably won't get much time to write about what I've been up to this weekend. So here's a little taster. My first robot programmed in Haskell:



It's a pretty simple device that's programmed to follow objects that come into the view of the three IR distance sensors mounted on the front of it. (The object in this case being a skull on the end of a stick, must have been left over from Hallowe'en...) But the point, at this stage, wasn't to do something sophisticated with the robot, but to see if I could actually program such a device in a high level language that wasn't C/C++ or assembler.

For hardware I used the Make controller and a Rigel robot base. I also use a bunch of Sharp IR sensors.

So for software I started with Matthew Naylor's module for Lava described in The Monad Reader Issue 7. I didn't have to do much with this code at all so it's almost embarassing to claim I did anything. The Make controller is ARM based and programmable in C. In fact it can run FreeRTOS and I simply modified Matthew's code to spit out a C function that is spawned as a task in FreeRTOS. It was about 30 minutes work to make the changes. With Matthew's code modified I now have a Haskell DSL that is translated to C for compilation for the robot. Strictly speaking there isn't actually a Haskell runtime on the robot. Think of the Haskell as a metaprogram that spits out C.

The actual Haskell DSL code looks something like this:

> main =  writeMakeController "rigel" f    where
>         f inp = Output {
>            servoSpeed0   = 545+leftSpeed,                       
>            servoSpeed1   = 550-rightSpeed,
>            servoSpeed2   = 550+leftSpeed,
>            servoSpeed3   = 540-rightSpeed
>          }
>             where
>                 rightSpeed = larger (threshold 300 leftSensor) (threshold 300 centreSensor)
>                 leftSpeed = larger (threshold 300 rightSensor) (threshold 300 centreSensor)
>                 centreSensor = smooth $ adc0 inp
>                 leftSensor = smooth $ adc5 inp
>                 rightSensor = smooth $ adc4 inp
> smooth x = (x + delay 0 x + delay 0 (delay 0 x) + delay 0 (delay 0 (delay 0 x)))/4
> threshold m x = x />=/ m ? (256,0)
> larger x y = x />=/y ? (x,y)

And the video above is the resulting behaviour.

Note, this is the first thing I got working. I make no great claims for the quality or ingenuity of the code. Most of it should be self-explanatory but note the function smooth which averages over the last inputs. This doesn't make sense when considered as an ordinary function because in a pure functional language you can't add the argument of a function to the value that you passed in last time you called the function. But these values are all streams and in principle this makes the DSL extremely powerful. Anyway, I'll write it up better at a future point and I'll see what I can do about making the code available.

Homeland Security Threat Level Monad

I've previously introduced the idea of a monad being a kind of one way wrapper. The Monad interface allows you to wrap things, but not unwrap them again. And it allows you to apply any function you like to a wrapped object, as long as it stays wrapped. I also indicated how this could be used to model a notion of taint. If you have some data from a tainted source - for example it's unreliable, then the trivial monad can do a nice job of keeping track of all other data that has been tainted by contact with this data. What I also want to show is that a wide variety of other monads can be interpreted in terms of taint. What this means is that despite its simplicity, the trivial monad is always worth keeping in the back of your mind because it is a protoype for many other monads.

You could imagine using this taint idea as a security model. Any data of type W a is considered secret and a threat to public safety if it is revealed. By sticking to the Monad interface we force the compiler to make sure that any data depending on secret data is itself marked secret. This makes it easy to keep secret data sequestered - at least as long as people only use the Monad interface to W.

But suppose this isn't good enough and we need to track not only which data is secret, but also the threat level to public safety should a secret be revealed. With this in mind, define the Homeland Security Threat Level Monoid:

> data Threat = LOW | GUARDED | ELEVATED | HIGH | SEVERE
>          deriving (Eq,Ord,Enum,Show)

It's a monoid because it has an associative binary operator max and has an identity LOW. For example max LOW ELEVATED = ELEVATED.

We want to tag some of our data with a threat level. In addition, we want to ensure that if we combine two tagged pieces of data, x and y, the combination must be tagged with the worse of the two threat levels. That's why we have the monoid structure that allows us to combine two threats.

So define

> data T a = T Threat a deriving Show

But to use this it appears we have to thread use of the max function all the way through our code. It'll be a nightmare to implement. But there's a solution. What we want is a monad!

Consider LOW to be the default threat level of all data and that tells us how to define return.

Now suppose we have a function f :: a -> T b and some data x :: T a. We want to be able to apply f to x using the monad interface: x >>= f. The threat level of the result should be the max of the two possible sources of threat: x itself or the result of computing f. With that in mind:

> instance Monad T where
>     return x = T LOW x
>     (T t x) >>= f = let T t' x' = f x in T (max t t') x' 

Notice how we've now shown how to implement threat combination in a way that can be used by any code that uses the monad interface. So now we can take last week's definition of h and use it:

> g x y = y >>= (return . (+x))
> h x y = x >>= (\x -> g x y)

Try evaluating h (T ELEVATED 3) (T GUARDED 6). It just works, even though g and h were originally written for the trivial monad! So it's trivial to make code that's written monadically work with our new threat level system. This should give some idea of the power of monads. We've completely changed the sematics of our monadic addition operator, h, without lifting a finger to edit it.

But notice that I didn't need to talk about threats at all. I've shown how to tag data with other data that's implicitly carried along with it. All we needed was a rule for combining pairs of tags. Our tags were of type Threat. Instead we could have used strings and instead of max we could have used ++. If we had, then we'd have a way of implicitly carrying around logs with our computations. In fact, all I have done here is implement a special case of the Writer monad. But I still think that thinking in terms of data that is quarantined in some way helps with the intuition.

I hope that by starting with the trivial monad last week, and upgrading in a natural way to the Homeland Security Threat Level monad, I've given an easy way to think about monads. Our ordinary intuitions about taint, for example that if a clean object touches a tainted one it becomes tainted itself, carry over to monads. And the Writer monad isn't the only one that can be thought of in these terms. I hope at some point soon to look at the Maybe and list monads in this way too.

The Trivial Monad

Update: I've added some solutions to exercises below.

Haskell monads form a type class. And Haskell type classes are essentially interfaces shared by some types. So the Monad type class is a common API for talking to a bunch of different types. So the question is this: what's so special about this API?

One way to grasp an API is to see concrete examples. Usually the monad API is introduced through containers or IO. But I think that the prototypical monad through which all others can be understood is much simpler - it's the trivial monad.

APIs often capture design patterns, and the design pattern here is a one-way wrapper. We can wrap things but we can't unwrap them. We still want to be able to whatever we like to the wrapped data, but anything that depends on anything wrapped should itself be wrapped.

Without further ado, here's a wrapper type:

data W a = W a deriving Show

(The Show bit is just so we can play with these things interactively.)

Note how it doesn't add anything except some wrapping. And the first thing we need is to be able to wrap anything we like. We can do that easily with this function:

return :: a -> W a
return x = W x

(We could have just used W instead of return. But I'm heading towards a common API that will work with other types too, so it obviously can't be called W.)

And now we need one more thing - a way to manipulate wrapped data leaving it wrapped. There's an obvious idea. Given any function a -> b we write a function that converts it to a function W a -> W b. This is guaranteed to keep things under wraps. So here goes

fmap :: (a -> b) -> (W a -> W b)
fmap f (W x) = W (f x)

.

It seems that we've finished. For example we can wrap a number and then increment it:

a = W 1
b = fmap (+1) a

But here's something we can't do. Define a function f like this:

f :: Int -> W Int
f x = W (x+1)

It increments and returns a wrapped result. Now suppose we want to apply the underlying operation here twice, ie. we'd like to increment a number twice and return a wrapped result. We can't simply apply f twice because (1) it'll doubly wrap our result and (2) it'll attempt to add 1 to a wrapped result, something that doesn't make sense. fmap doesn't do what we want either. Try it in an interactive Haskell session. It seems we need a way to apply f, unwrap the result and apply f again. But we've already said that we don't want to allow people to unwrap these things. We we need to provide a higher order function that does the unwrapping and application for us. As long as this function always gives us back something that is wrapped, our end users will never be able to unwrap anything. Here's an idea for such a function:

bind :: (a -> W b) -> (W a -> W b)
bind f (W x) = f x

Notice how it's very similar to fmap but is even simpler. So now we can try doubly, or triply, incrementing

c = bind f (f 1)
d = bind f (bind f (f 1))

Notice how bind is more general than fmap. In fact, fmap f = bind (return . f).

And that's it. Using return and bind we have achieved our goal of wrapping objects and freely manipulating wrapped objects while keeping them wrapped. What's more, we can chain functions that wrap without getting bogged down in multiple layers of wrapping. And that, really, sums up what a Haskell monad is all about.

So here are a couple of exercises:

(1) define a function g :: Int -> W Int -> W Int so that g x (W y) = W (x+y). Obviously that definition won't do - the left hand side has a W y pattern so it's actually unwrapping. Rewrite this function so that the only unwrapping that happens is carried out by bind.

(2) define a function h :: W Int -> W Int -> W Int so that h (W x) (W y) = W (x+y). Again, no unwrapping.

I'm hoping that after you've done these exercises you'll see how you can still work freely with data even though it's wrapped.

In Haskell, it is more usual to use the operator >>= instead of bind where bind f x = x >>= f.

So the last question is this: why would you ever wrap data like this? In practice people tend not to use the trivial monad very much. Nonetheless, you can see how it might be used to represent tainted data. Wrapped data is considered tainted. Our API never lets us forget when data is tainted and yet it still allows us to do what we like with it. Any time we try to do anything with tainted data the result is also tainted, exactly as we might expect. What I find interesting is that almost every monad, including IO, lists and even probability, can be thought of quite naturally as variations on taint. I hope to say more about this in the near future.

Anyway, code above fails to compile because of some namespace clashes. Here's a complete definition of W that really works. Note also that fmap, which we showed we don't really need, allows us to make W an instance of Functor too:

> data W x = W x deriving Show

> instance Functor W where
>    fmap f (W x) = W (f x)

> instance Monad W where
>    return x = W x
>    W x >>= f = f x

Exercise 3: Prove the three monad laws for W. This should be almost trivial.

Exercise 4: We can't completely unwrap things using the monad API. But we can unwrap one layer from things that are wrapped twice. So here's a nice puzzle: define a function join :: W (W a) -> W a using the Monad API and no explicit unwrapping.

In conclusion, most monad tutorials show what monads are, and how to use them. I hope I've given some additional insight into why the Monad interface consists precisely of the two functions it does. That insight should become clearer when I say more about taint.

PS I've been a bit busy lately with planning for kitchen remodelling, lots of work and a new cat(!). I probably won't be posting very frequently for a few months.

---

Solutions to Exercises

Some people have found the exercises tricky so I'm posting some solutions. I think it's a good sign that they're tricky, it means that there is some non-trivial stuff to be learnt from playing with this monad. It ought to figure more strongly in tutorials. In a sense it allows people to learn about monadicity in its purest form.

Anyway, enough of that.

Exercise 1: We want g :: Int -> W Int -> W Int so that g x (W y) = W (x+y). Start the definition like this

g x y = ...

We want to apply (+x) to y. But y is wrapped so we must use >>= to apply it. But (+1) doesn't have the right type to be applied with >>=, it's an Int -> Int not an Int -> W Int. So we must compose with return and we get

g x y = y >>= (return . (+x))

Exercise 2: We want h :: W Int -> W Int -> W Int so that h (W x) (W y) = W (x+y). But g is already fairly close. So again write

h x y = ...

The difference from g is that x is wrapped. So we want to apply \x -> g x y to our wrapped value. This function is already of type Int -> W Int. So we can just write

h x y = x >>= (\x -> g x y)

or just h x y = x >>= flip g y.

Exercise 3: I'll just do the last one: (m >>= f) >>= g = m >>= (\x -> f x >>= g)


m is of the form W x for some x:

(W x >>= f) >>= g 
= (f x) >>= g
 
W x >>= (\x -> f x >>= g)
= (\x -> f x >>= g) x
= f x >>= g

So the LHS equals the RHS.

Incidentally, you can get a long way with these problems by not thinking about what's going on! Instead, just write out the type signatures of all of the functions involved and try to stick them together like a jigsaw puzzle so that the final result has the right type. At each stage the arguments of each function must match the signature of the function so it really is like fitting the shapes of jigsaw pieces together. This isn't a foolproof strategy, but it often works.

Anyway, feel free to ask questions if you're still stuck. Some of your questions may be answered here.

The curious rotational memory of the Electron, Part 2

A spinning top keeps going because it has angular momentum, and the angular momentum arises because the top is made up of mass that is distributed around the axis of rotation. The angular momentum is represented by a vector (well, a pseudovector really) directed along the axis of spin (at least when it's spinning around the top's axis). Similarly electrons have angular momentum. It's tempting to think of electrons as like microscopic spinning tops but the angular momentum of electrons isn't due to an extended mass rotating about some axis. In fact, as far as we know, electrons have no extension. Instead their angular momentum is a peculiar kind of intrinsic angular momentum that has a lot in common with ordinary angular momentum - eg. it's conserved and it's interconvertible with ordinary angular momentum. Because of it's not being associated with a mass distribution it's given the special name of spin, but as long as you bear in mind these provisos it's probably fair to use at least some intuitions about spinning tops when thinking about electrons.

But it's a little more complex than this. We also have quantum mechanics to contend with. The spin of an electron is a vector. But we find that when we measure one of the components of this vector this value is quantised and can only take values +hbar/2 and -hbar/2, where hbar is Planck's constant. We choose units where h-bar is 1 so the z-component of the spin is always measured to be +1/2 or -1/2. If we write these two states as |+> and |-> then because we are dealing with quantum mechanics, the z-component of the spin can be represented by the linear combination a|+>+b|->. This corresponds to a state in which there is a probability |a|² of measuring +1/2 and a probability |b|² of measuring -1/2. This is what might have been written as a.*return (1/2)+b.*return (-1/2) in my earlier Haskell code. But that's just one component of the spin. What about the x- and y-components? Amazingly the state a|+>+b|-> tells us everything we can possibly know about the spin of an electron and we'll call it a spin state.

Suppose we have an electron in the state ψ = a|+>+b|->. What happens if we measure the y-component of its spin? One way to answer that question is to rotate the electron through π/2 so that its x-axis is rotated to align with the z-axis and then measure the z-component of its spin. In order to to that we need to know how to rotate spin states. The rule for rotation through θ about the x-axis is this (in a suitable coordinate frame):

|+> → cos(θ/2)|+>-sin(θ/2)|->
|-> → sin(θ/2)|+>+cos(θ/2)|->

Note how choosing θ=0 gives the identity, as expected. Note also that θ=π maps a|+>+b|-> to b|+>-a|-> so that the probabilities of measuring +1/2 and -1/2 are simply swapped, exactly what you'd expect for turning a state upside down. But there's something else that you should notice - there's an ambiguity. A rotation through 2π should give the same as a rotation through 0 and yet setting θ=2π in that transformation maps a state ψ to -ψ. Now |a|² = |-a|² so the probability of observing spin up or spin down is unaffected. But as I've been showing over previous posts, flipping a sign in a state can make a big difference as soon as you start performing interference experiments. The same goes for any angle: if I rotate through π should I use θ=π or θ = 3π? So can the transformation I've given make sense?

The transformation does make sense if you consider that in any physical process that rotates an electron the transformation will evolve continuously over time. Electrons don't just instantly rotate. In other words, if a rotation is applied to an electron then it will follow a path in SO(3), not just be an instantaneous application of an element of SO(3). And that allows us to resolve the ambiguity: the rotations of electrons are described by the double cover of SO(3) known as SU(2). So a rotation through 360 degrees doesn't return you to the identity although a 720 degree rotation does. The transformation I gave above is completely unambiguous if you continuously rotate an electron around the x-axis tracking a continuous value of θ, after all, the double cover is basically just the set of continuous paths from the identitiy in SO(3) (with homotopic paths considered equivalent).

And that's the bizarre fact: electron rotations aren't described by SO(3), they're described by SU(2). In particular, rotating an electron through 360 degrees does not return it to its original state, but a rotation through 720 degrees does! In a sense, like Dirac's belt, electrons can remember something about the path they took to get where they are, in particular they remember how many twists there were in the path.

What does this mean experimentally? the first thing to note is that this is true not just for electrons but any spin-1/2 fermion. This included protons and neutrons. The stuff I've been talking about manifests itself in a number of ways. In particular, the spin of a particle affects how a magnetic field acts on it. For example, spin-up and spin-down particles can be separated into distinct beams using Stern-Gerlach apparatus. Also, the spin of particles precesses in a magnetic field and this is used on a regular basis in NMR. These two facts allow us to easily manipulate and measure the spin of fermions. In other words, the fact that fermions remember how many twists there are in their rotations isn't just some esoteric nonsense, it's now engineering and the theory is tested repeatedly all over the world.

Every familiar object is invariant under rotations through 360 degrees. So the fact that electrons need to be rotated through 720 degrees to return them to their original state seems like one of the most bizarre facts about the universe I know of. And yet many books that introduce spin just slip in this fact in a routine way as if it were no different to any other.

The fact that the biggest connected cover of SO(3) is the double cover puts a big constraint on the kinds of weird effects like this can happen. We can have a 360 degree rotation multiply by -1, but not by i, because a 720 degree rotation absolutely has to return us to where we started from. But suppose the universe were 2-dimensional. If you remember what I said about SO(2) you may notice that no such constraints apply because SO(2) has an infinite cover. There is a group in which all of the rotations through 360n degrees are distinct for distinct n. This means that a physical system could have its state multiplied by any factor (of modulus 1) when rotated through 360 degrees. Particle that behave this way are called anyons. But we live in a 3D universe so we don't expect any fundamental particles to have this property. However, in quantum mechanics any kind of 'excitation' of a physical system is quantised and can be thought of as a type of particle. These are known as quasiparticles. For example, just as light is made of photons, sound is also quantised as phonons. In the right kind of solid state medium, especially those that arise from some kind of 2D lattice, it seems quite plausible that anyons might arise. This gives rise to the so called fractional quantum hall effect. Anyons might one day play an important role in quantum computing via topological quantum computation.

The Curious Rotational Memory of the Electron, Part 1

There's a curious and bizarre fact about the universe that is introduced in physics courses without anyone stopping to point out just how curious and bizarre it is. I think it deserves to be more widely known to non-physicists. It's a little esoteric, but not so esoteric it can't be grounded in ordinary everyday concepts. And it isn't something hypothetical but something that is fundamental to modern physics and testable in the lab. It is also the reason why I chose to use theta/2 instead of theta in my quantum computing code.

This post will be in two parts. This part will mainly be about geometry and topology, and in the second part I'll talk about the physics.

In order to explain the idea I want to introduce the concept of a universal cover. If you look at the definition of this notion anywhere it can seem pretty daunting to non-mathematicians, but in actual fact the key idea is so simple that even a child could understand it. I can say that without fear of exaggeration because universal covers used to play a role in a childhood fantasy of mine. (OK, I was a weird kid, but not that weird.)

My crazy idea as a kid was this: if there were two ways to get from A to B, maybe the B you reached by taking either of the two routes were actually slightly different, even though they looked the same. I was actually more concerned with journeys from A to A. Suppose the kitchen has a door at each end and that you start in the bedroom, go down to the kitchen at one end, come out the other end, and return by a different path to the bedroom. How do you know this is the same bedroom? (I also admit I read a lot of science fiction as a kid.)



Suppose instead you took a simpler path: you left the bedroom and went to another room and then returned by the path you went. If you had taken a piece of string and tied one end to the bedroom and walked this path playing out the string behind you as you went then when you returned to your bedroom you would have both ends of the string visible. You could then reel in the string from your journey, without letting go of the ends, and see that indeed, both ends were in fact the ends of the same string and so you had retuned to the start of your journey. But if you'd taken the round trip through the double-doored kitchen you wouldn't be able to pull the string tight and so you could no longer be sure.

In what follows we consider two paths from A to B equivalent in this way: mark the two paths with strings going from A to B. If you can pull one string to line up exactly with the other (you're allowed to stretch) without breaking the string or moving the endpoints then they are considered equivalent. (The correct technical term is homotopic.) So if a path that appears to go from A to A is homotopic to the trivial zero length path, then you know that the path really is from A back to A. But if you can't shrink the path down to zero length you're faced with the possibility that the string is extended from A to a similar looking place, A'. (Of course this implies a considerable conspiracy with someone tying a distinct identical string to an identical bedroom, but if you can't see both ends, alien beings from a higher dimension could do anything...)

The best case scenario is that all of the bedrooms are the same - a conventional house. Call this house X. But what's the worst case scenario? How many different copies of the bedroom might there be? Well let's call the starting bedroom 0. If we walk anti-clockwise round the house, going through the kitchen, we end up back at a bedroom. In the worst case scenario this will be another bedroom. Let's call it 1. Walk round again and we get to another copy of the bedroom, call it 2 and so on. So we have an infinite sequence of bedrooms, (0,1,2,…). But we could have walked clockwise. In that case we'd get a sequence (-1,-2,…). So now think about the set of all points in this infinite house. Each point can be described as a pair (x,n) where x is a point in the house X and n labels which of the infinite number of paths you could have taken to get there. (Remember, if the two paths are homotopic we consider them to be the same path.)

If using a house as an example seems too far-fetched think of a mult-storey car park. It's easy to walk what you think is a closed loop but find that you've in fact walked to almost identical looking, but different, floor.

(Incidentally, this kind of infinite house is the kind of thing that can be created in a Portal engine such as in the game Portal.)

Let's think about this more generally. Given a topological space X we can form another space X' from it by (1) picking a point p in it and (2) forming the set of all pairs (x,n) where x is a point in X and n is a path (up to homotopy) from p to x. X' is known as the universal cover of X.

Enough with houses, let's try a more purely geometric example. Consider the space of rotations in 2D, otherwise known as SO(2). A rotation can simply be described by an angle with the proviso that a rotation through an angle θ is the same as a rotation through θ+2nπ where n is an integer. The space of rotations is essentially the same thing as the circle (by circle I mean the 1-dimensional boundary of a 2-dimensional disk). Looping once round the circle brings you back to where you started. So how can we visualise the universal cover of SO(2)? There's a nice trick for this. Imagine a disk that's free to rotate around one axis. SO(2) describes the rotations of this disk. But now imagine a belt, fixed at one end, connected at the other end to the centre of this disk like this:



The rule is that we're allowed to twist the belt as much as we like as long as the centre line of the belt is along the axis. Notice how if we rotate the disk through an angle of 2nπ the disk returns to its starting point, but the belt 'remembers' how many twists we took to get there. Notice another thing: if we travel along the belt from one end to the other, the orientation of the belt at each point gives a 2D rotation. In other words, the belt itself represents a path in SO(2) with one one representing a fixed point in SO(2) and the other end matching the current orientation of the disk. So if the orientations of the disk represent elements of SO(2), the possible twistings of the belt represent elements in the universal cover of SO(2). And what the twists do is allow us to differentiate between rotations through θ and rotations through θ+2nπ. In other words, the universal cover is just the real line, but not folded down so that θ is the same as θ+2nπ.

The universal cover is the worst case scenario where all non-homotopic paths are considered to represent distinct points. But it's also possible to construct an in-between scenario where some paths that aren't homotopic to each other are still considered to lead to the same place. For example, take the real line again. If we identify 0 and 2π we get SO(2). If we consider 0 and 2nπ to be distinct for n≠0 we get the universal cover. But we pick an N≠1 and consider 2Nπ to be the same as 0, but that 2nπ is distinct from 0 for 0<n<N. What we now end up with is an N-fold cover of SO(2). Unfortunately I don't know a way construct a physical example of a general N-fold cover using belts and discs. (At least for N>2. For N=2 you should easily be able to come up with a mechanism based on what's below.)

With the examples so far, the universal cover of a space has always had an infinite number of copies of the original space. That isn't always the case and now we'll meet an example where the universal cover is a double cover: the real projective plane, RP². The easiest way to visualise this is to imagine the game of Asteroids played on a circular TV screen. In the usual Asteroids, when you go off one side of the screen you reappear on the opposite side. In this version it's much the same, points on the edge that are diametrically opposite are considered to be the same points. Such a universe is a real projective plane. Let P be the point at the centre of the screen. Suppose we attach one end of our string to P and fly off towards the edge, reappearing on the other side. Eventually we get back to P again. Can we shrink this path down to a point? Unfortunately, any path that crosses the edge of our screen will always cross the edge no matter how much we try to deform it. If we move the point at which it crosses one edge of the screen, the point at the other side always remains diametrically opposite.




This means you can never bring these points together in order to collapse down a loop. On the other hand, suppose we loop twice round the edge. The following diagram shows how this path can be shrunk down to a loop:



So given two distinct points in RP² there are two inequivalent paths between them, and so the universal cover of RP² is made from two copies of the space and is in fact a double cover.

Now we'll try a 3D space - the space of 3D rotations, otherwise known as SO(3). This space is 3-dimensional because it takes 3 numbers to define a rotation, for example we can specify a rotation using a vector whose direction is the axis of rotation and whose length is the angle of rotation. Now all possible rotations are rotations through an angle 0≤θ≤π for some exis. So every rotation can be represented as a point in the ball of radius π. But a rotation around the axis a through an angle of π is also the same as a rotation around the axis -a through π. In other words, SO(3) is the closed 3-dimensional ball with diametrically opposite points identified. This is very similar to the situation with RP². In fact, almost exactly the same argument shows that the universal cover of SO(3) is the double cover. (For an alternative wording, see this posting.) There is also another way of seeing this that is analogous to the disk with belt example above. This time we have a sphere instead of a disk. The belt is fixed at one end and the other end is attached to the sphere. We are free to rotate the sphere however we want, and we're no longer confined to having the belt along one axis. It is an amazing but easily demonstrable fact that if you rotate the sphere through 2π, you cannot untwist the belt, but if you rotate it through 4π you can! This is known as the Dirac Belt Trick, and if you don't believe it, watch the Java applet written by science fiction author Greg Egan.

If you've done much computer graphics you'll have already met the double cover of SO(3). Computer graphics people sometimes like to represent rotations using unit quaternions. The unit quaternions form a 3-dimensional sphere but there are two such quaternions for each rotation. So the unit quaternions actually represent this double cover, and SO(3) is a 3-dimensional sphere with antipodean points identified.

And that's the geometry dealt with. In the next installment I'll say what this has to do with electrons and my quantum computing code.

---

Footnote

Balls and Spheres. The closed n-dimensional ball is the set of points (x₁,...,x_n) such that x₁²+...+x_n²≤1. The (n-1)-dimensional sphere is the set of points (x₁,...,x_n) such that x₁²+...+x_n²=1.

The Angel Problem has been Solved! (Maybe)

Forget computing the Kazhdan-Lusztig polynomials for E8, a far more important theorem has apparently been proved. I won't state the answer here as it's fun to think about the problem for yourself. But a number of proposed proofs of solutions to the Angel problem have been informally published, at least electronically.

Firstly, what is the problem? Wikipdia has a succinct write up. But here's a quick summary anyway:

The game of Angels and Devils, invented by John Conway (the same Conway as in the Game of Life, the Monstrous Moonshine conjectures, the proof that subatomic particles have free will and so on), is played on an infinite chessboard. Fix some constant k. One player plays the Angel. The other plays the Devil. On the Devil's turn, he may eat one square on the chessboard. (The devil doesn't actually have a piece, he just eats whatever square he wants.). On the Angel's turn she jumps to a square (distinct from the current one) up to k king's moves away. The Devil wins if the Angel has no legal move. The Angel wins if she has a strategy that guarantees always to have a legal move. For each k, who wins? Until recently it was unknown if the Angel could win even for arbitrarily large k. This may seem like a surprise, after all, if the angel can move 1,000 king's moves, surely she can easily evade the jaws of the devil who can only chomp one square at a time. But surprisingly, the obvious Angel strategies that come to mind can be countered.

What's nice about this problem is that it's so easy to state, and yet has stumped some of the smartest people who like to think about these problems.

Two papers with possible solutions are Máthé's and Kloster's. I don't think either of them have formally been published yet so I don't know if the peer review is complete.

The Shor Quantum Error Correcting Code (and a Monad for Heat)

Quantum computers promise great things, but that promise will only come to be if quantum computers can be made reliable. In order to achieve that reliability we need qubits that can maintain their state for extended periods of time. Unfortunately, just about any type of interaction between a qubit and its environment can render the qubit useless. But there is a solution - a quantum error correcting code. This is very similar to classical error correcting code in which a number of bits are stored, in a distributed fashion, over a larger set of bits, in such a way that if some of the bits are corrupted, the original state can be recovered. In the following will be some code to simulate the implementation of one of these codes: the nine bit Shor code.

Firstly, let me include the quantum monad that I've developed over earlier weeks:

> import Prelude hiding (flip)
> import Maybe
> import List
> import Data.Map (toList,fromListWith)
> import Complex
> import Data.Bits hiding (rotate)

> infixl 1 ->-, -><, =>=
> infixl 7 .*

Remember that W b a is the type of elements of a vector space over b with a basis whose elements are labelled by elements of a. I'm using the letter W as I'm thinking of elements of b as being weights.

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

It's useful to be able to apply transforms to all of the weights:

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

W is an instance of Num, Functor and Monad. Vector spaces form monads because you can specify a function between vectors spaces in terms of how it maps basis elements in the first vector space. The monad bind function lifts the function on basis elements to the entire vector space.

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

(.*) gives us scalar multiplication:

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

It'd be cool if the bind function for W could simplfy our vectors by collecting together multiples of the same basis elements. Unfortunately this is tricky to do for reasons described in Eric Kidd's recent Set monad article. Instead I'll use an explicit function, collect, that collects our terms together. I've also used a crude method for throwing out terms that probably ought to be zero. A better solution might be to work only with rational numbers, or exact algebraic numbers.

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

And here is the quantum monad:

> type Q a = W (Complex Double) a

There are a number of ways to apply functions to vector spaces. I'll be writing functions on the right so that I can apply operations to a quantum state reading from left to right. So here's a bunch of combinators to implement this. First there's ordinary function application:

> (->-) :: a -> (a -> b) -> b
> x ->- f = f x

Then there are functions that operate on the vector space by simply relabelling the basis:

> (=>=) :: Q a -> (a -> b) -> Q b
> x =>= f = fmap f x

Then there are functions that map from Q a to Q b but are defined by how they act on the basis elements of Q a:

> (-><) :: Q a -> (a -> Q b) -> Q b
> g ->< f = g >>= f

As I mentioned last post on entanglement, we can form states for joint independent systems using the tensor product:

> (<*>) :: Q a -> Q b -> Q (a,b)
> x <*> y = do
>   a <- x
>   b <- y
>   return (a,b)

The tensor product is a lot like a multiplication table. The coefficient of (a,b) in x <*> y is the product of the coefficients of a and b in x and y respectively. If two systems are unentangled then you can simulate the physics of each part separately. So it'd be nice to have a way to test whether or not a given state is unentangled, in other words, to test whether a state is a tensor product. As I'm not developing a serious application here, here's a quick and dirty algorithm that more or less does what anyone would do when faced with the task of recovering the row and column headings from a multiplcation table:

> norm x = sqrt $ sum $ map (\(_,a) -> magnitude a^2) $ runW $ collect x

> untensor :: (Ord a,Ord b) => Q (a,b) -> (Q a,Q b)
> untensor (W a) =
>     let u = nub $ sort $ map (fst . fst) a
>         v = nub $ sort $ map (snd . fst) a
>         p = W $ map (\x -> (x,fromJust (lookup (x,head v) a))) u
>         q = W $ map (\x -> (x,fromJust (lookup (head u,x) a))) v
>         in ((1/norm p :+ 0) .* p,(1/norm q :+ 0) .* q)

untensor recovers x and y from x <*> y (up to multiplication by a complex number with unit norm) but it gives garbage if the argument isn't a tensor product. So we need to test whether or not the original vector is the tensor product of the inputs. So here's a test of how close the original is to the product of the putative factors:

> independence x = collect (x-uncurry (<*>) (untensor x))

We'll be doing lots of stuff with tensor products so it's useful to have some tools to apply functions to one or other factor:

> mapFst f (a,b) = (f a,b)
> mapSnd f (a,b) = (a,f b)
> mapFstM f (a,b) = do
>     a' <- a ->- f
>     return (a',b)
> mapSndM f (a,b) = do
>     b' <- b ->- f
>     return (a,b')

If a bit is a Bool then by rights, a qubit should be a Q Bool. But I think 0 and 1 look nicer so I'm using Q Integer instead. So here's the fundamentally quantum rotation operator:

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

Although rotate theta is a map Integer -> Q Integer, using ->< it becomes a map Q Integer to Q Integer. Note that the Double argument isn't float data inside a quantum device, it's a hard-coded feature of the hardware itself. For example, if we choose theta to be pi/2 and -pi/2 we get the standard square root of not operation and its inverse:

> hnot = rotate (pi/2)
> unhnot = rotate (-pi/2)

(You may be wondering about the use of theta/2 in my definition of rotate. There's an interesting story behind that involving what I think is one of the weirdest features of physics. For another day.)

I'm also goint to manipulate lots of lists of qubits so here are some tools for that:

> mapAt i f x = let (h,t:ts) = splitAt i x in h ++ (f t : ts)

> mapAtM :: Monad m => Int -> (a -> m a) -> [a] -> m [a]
> mapAtM i f x = let (h,t:ts) = splitAt i x in do
>     y <- f t
>     return $ h ++ (y : ts)

> flip a = 1-a
> flipAt i = mapAt i flip

So here's a statement of the problem. We have a qubit, ie. something of type Q Integer. But quantum systems are very delicate and our qubit may wander into a different state. In particular, after a time t it will be as if rotate theta has been applied to it, except that we can't control theta. We need to find a way to recover the original qubit. In fact, it's worse than that: our qubit is in some environment which we can consider to be a list of qubits. So our state is given by an instance of Q ([Integer],Integer) and the interaction might not simply be a rotation of the qubit but an operation that entangles the qubit with the environment. This might happen if the qubit isn't perfectly isolated.

Let's forget about the full generality here and consider just one type of data corruption - a qubit that might or might not have been flipped. How can we make a device robust against such occurences? Well this is practically identical to the problem of making classical data robust against noise. We can simply use a 3 bit repetition code to store our qubits as three qubits. So here's an implementation of a function to repair a codeword from the classical 3 bit repetition code:

> repairRepetition' [b,c,d] = (if a<3 then mapAt a flip else id) [b,c,d] where
>     a = case (b,c,d) of
>         (0,0,0) -> 3
>         (1,1,1) -> 3
>         (0,0,1) -> 2
>         (1,1,0) -> 2
>         (0,1,0) -> 1
>         (1,0,1) -> 1
>         (1,0,0) -> 0
>         (0,1,1) -> 0

But there are a couple of problems with this. Firstly this is a many-to-one map because [0,0,0] and [0,0,1], say, both get mapped to [0,0,0]. In quantum mechanics, all of our maps need to be unitary, and in particular, invertible. A many-to-one map isn't invertible. Secondly, we're creating a register a, and again this is not an invertible operation. These two problems will actually cancel each other out because the apparent destruction of information that comes about from our many-to-one map can be fixed by allowing a to carry that information away. And we can ensure reversibility by using a supplied register rather than creating a new one:

> repairRepetition'' (a,[b,c,d]) = (a,(if a<3 then mapAt a flip else id) [b,c,d]) where
>     a = case (b,c,d) of
>         (0,0,0) -> a `xor` 3
>         (1,1,1) -> a `xor` 3
>         (0,0,1) -> a `xor` 2
>         (1,1,0) -> a `xor` 2
>         (0,1,0) -> a `xor` 1
>         (1,0,1) -> a `xor` 1
>         (1,0,0) -> a `xor` 0
>         (0,1,1) -> a `xor` 0

It's easy to write an inverse for this function so it's a valid unitary transform. But it's actually a bit tedious to write out all those xors. So we can cheat a bit. Let's just assume there's an infinite supply of input zeroes and use them at our leisure. We no longer need to write the xors, we just note that every time we assign a value to a new register the hardware engineer building your device installs some means of import zero qubits. The a register in the above function will never be touched again. But we can't simply throw it away. It's a kind of side effect and we'll call these bits <em>ancillae</em>. We'll have lots of these side effects that we need to collect up, so here's the final quantum repetition codeword repair code:

> repairRepetition (a,[b,c,d]) = (i:a,(if i<3 then mapAt i flip else id) [b,c,d]) where
>     i = case (b,c,d) of
>         (0,0,0) -> 3
>         (1,1,1) -> 3
>         (0,0,1) -> 2
>         (1,1,0) -> 2
>         (0,1,0) -> 1
>         (1,0,1) -> 1
>         (1,0,0) -> 0
>         (0,1,1) -> 0

There's one last thing we could do to this code: the pairs being handed in and out of this function are essentially just instances of the Writer monad. So we could just squirrel away our ancillae inside a monad. I won't do that here because I think it's good to see everything explicitly when figuring this stuff out for the first time. But think about the bits that we'd be writing out. They correspond to information, or entropy, that we'd ultimately like to pump out of our device as exhaust. In other words they are heat. So this monad would, in effect, be a heat monad. I found this amusing given what I wrote in reason number eight in my list from last year.

So far so good. But that only corrects one type of error: a bit flip. But note that it also protects superpositions from bit flips. So let's take a quantum state built from repetition codewords placed in an empty environment:

> state1 = return [] <*> (return [0,0,0] + return [1,1,1])

Now corrupt the last bit:

> state1' = state1 =>= mapSnd (mapAt 2 flip)

And repair it:

> state1'' = state1' =>= repairRepetition

Although we have repaired the state of our system, the state of the environment has changed. But that's fine because this state is still an unentangles tensor product of the environment and our system. You can see this just by looking at the state, but for peace of mind use independence state1'' to see that the deviation from being a tensor product is the zero vector.

But that's still only a narrow class of error that we can fix. Another type of error maps return 0 to return 0 but maps return 1 to -return 1. This is known as a phase flip. One way to view this as as a flip in a rotated basis:

> phaseFlip x = x ->- hnot ->- fmap flip ->< unhnot
> pFlipAtM i = mapAtM i phaseFlip

There's a nice strategy to protect against phase flips: rotate our data into a basis where phase flips become ordinary flips, and use the repetition code. So we'll use these two states to represent 0 and 1:

> plus,minus :: Q [Integer]
> plus = return [1,1,1] ->< mapM unhnot 
> minus = return [0,0,0] ->< mapM unhnot 

So let's make a state that is a superposition of these states in some environment:

> state2 = return [] <*> (plus + minus) ->- collect

Corrupt the middle bit with a phase flip:

> state2' = state2 ->< mapSndM (mapAtM 1 phaseFlip)

And repair it:

> state2'' = state2 ->< mapSndM (mapM hnot) =>= repairRepetition ->< mapSndM (mapM unhnot) ->- collect

Comparing state2 and untensor state2'' we see that the repair is successful. Unfortunately, we can use this scheme to repair phase flips, but not plain flips. We need a scheme that can repair both of these. So now we turn to Shor's ingenious solution: use both of these schemes as the same time. So we'll use the scheme that protects against phase flips with the twist that instead of using three qubits we use three triplets of qubits. We can use the three triplets to fix phase flips and we can use the three bits within a triplet to fix flips. We'll represent 0 and 1 using shor 0 and shor 1 as defined here:

> tripletPlus = (1/sqrt 2) .* return [0,0,0] + (1/sqrt 2) .* return [1,1,1]
> tripletMinus = (1/sqrt 2) .* return [0,0,0] - (1/sqrt 2) .* return [1,1,1]

> shor 1 = do
>     a <- tripletPlus
>     b <- tripletPlus
>     c <- tripletPlus
>     return $ a ++ b ++ c

> shor 0 = do
>     a <- tripletMinus
>     b <- tripletMinus
>     c <- tripletMinus
>     return $ a ++ b ++ c

Repair is a little more complex before. First we use 6 ancilla bits to fix qubits within each triplet using the repetition code three times:

> ancilla (i,j,k) x = ((x!!i) `xor` (x!!j),(x!!j) `xor` (x!!k))

> fix (i,j,k) (e,x) = let (u,v) = ancilla (i,j,k) x
>                         b = case (u,v) of
>                             (0,0) -> x
>                             (0,1) -> flipAt k x
>                             (1,0) -> flipAt i x
>                             (1,1) -> flipAt j x
>                      in (u:v:e,b)

> flipRepair start = start
>             =>= fix (0,1,2)
>             =>= fix (3,4,5)
>             =>= fix (6,7,8)

We also use two ancilla bits to locate phase flips between the triplets. phaseFlipFix will be applied in a rotated basis inside phaseFlipRepair:

> phaseFlipFix :: ([Integer],[Integer]) -> Q ([Integer],[Integer])
> phaseFlipFix (e,y) = do
>     let u = (y!!0) `xor` (y!!1) `xor` (y!!2) `xor` (y!!3) `xor` (y!!4) `xor` (y!!5)
>     let v = (y!!3) `xor` (y!!4) `xor` (y!!5) `xor` (y!!6) `xor` (y!!7) `xor` (y!!8)
>     let c = case (u,v) of
>             (0,0) -> y
>             (0,1) -> y ->- flipAt 6
>             (1,0) -> y ->- flipAt 0
>             (1,1) -> y ->- flipAt 3
>     return (u:v:e,c)

> phaseFlipRepair x = x
>   ->< mapSndM (mapM hnot)
>   ->- collect
>   ->< phaseFlipFix
>   ->- collect
>   ->< mapSndM (mapM unhnot)
>   ->- collect

Now we can pull all of this together into a single process to repair both kinds of error:

> shorRepair x = x ->- flipRepair ->- collect ->- phaseFlipRepair ->- collect

So now all that's needed is to try corrupting some qubits to see if we can repair them.

> state3 = return [] <*> shor 1
> state3' = state3 ->< mapSndM (pFlipAtM 0)
> state3'' = state3' ->- shorRepair

Let's try a more complex superposition with a more complex corruption:

> state4 = return [] <*> (0.6.*shor 0 + 0.8.*shor 1) ->- collect
> state4' = state4 ->< mapSndM (pFlipAtM 5) =>= mapSnd (flipAt 5) ->- collect
> state4'' = state4' ->- shorRepair ->- collect

This is all well and good. But there are more ways to corrupt a qubit beyond simple flips. Let's mess up our state using two additional kinds of 'noise' that we haven't discussed. Firstly, instead of just flipping, we'll rotate one of the bits. Secondly, we'll additionally 'leak' that bit to the environment so as to cause decoherence. We've discussed solutions to neither of these problems, so let's see what happens.

> leak i x = do
>   (e,v) <- x
>   return ((v!!i):e,v)

> state5 = state4
> state5' = state5 ->- leak 4 ->< mapSndM (mapAtM 4 (rotate (pi/10))) ->- collect
> state5'' = state5' ->- shorRepair ->- collect

A minor miracle has taken place. The Shor code is robust against just about anything you can throw at a single qubit. In fact, any code that can survive flips and phase flips of a single qubit can survive any kind of corruption, including entanglement with the environment. (Despite having read the proof, I was still pretty surprised when the code actually worked!) And that gives researchers one hope that these methods might serve as ways to make quantum computing reliable.

I'll mention a couple more things in closing. Firstly, nowhere above have I mentioned wavefunction collapse, observations or superoperators. You don't need any of these things to define the Shor code and show that it can correct errors. (You probably do need some extra assumptions about quantum mechanics if you want to calculate probabilities of various outcomes in the event that the code fails - eg. due to multiple qubits being corrupted.) And secondly, despite it being my April 1st post last year, I still have the doubts I expressed here about quantum computing.

Oh...and expect those lines to take a few seconds to evaluate. When you simulate the Shor code you're performing computations in a 512-dimensional vector space!

Bibliography

1. Peter Shor's original paper has good discussion of the issues but I found the description of the correction algorithm a little laconic.
   
2. John Preskill's Lecture Notes were probably the most useful document.
   
3. A student term paper by Dheera Venkatraman was also pretty useful though I was confused for a while by what might be considered errors in it.
   

Update
Fixed a couple of minor typos. And remember this code should just run, without modifications, in GHCi, so you can play with it yourself.

Independence, entanglement and decoherence with the quantum monad

Before I can do anything more sophisticated with my code for probability theory and quantum mechanics I need to talk about a couple of related subjects - independence and entanglement. Again, I'm using the fact that formally, quantum mechanics is almost identical to probability theory, to kill two birds with one stone. So firstly, here's the probability/quantum code that I need from before:

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

> 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 Float a
> type Q a = W (Complex Float) a

Now consider this probability distribution:

> dice = sum [1/36 .* return (i,j) | i <- [1..6], j <- [1..6]]

This corresponds to the usual distribution you expect from rolling two dice. As we know, under normal circumstances, the outcome of one die is doesn't depend on the other - in other words they are independent. We can express this independence mathematically as follows. First make the distribution for one die:

> die = sum [1/6 .* return i | i <- [1..6]]

Now we can write code to make two independent die rolls:

> two_rolls = do
>   a <- die
>   b <- die
>   return (a,b)

We find that dice=two_rolls. We can make this more generic by defining a joint distribution as follows:

> tensor dist1 dist2 = do
>   a <- dist1
>   b <- dist2
>   return (a,b)

So tensor die die=dice. (This operation is actually a tensor product.) So now we can give a definition: if we have a distribution on (a,b), then the first and second components are said to be independent if the distribution can be written as tensor dist1 dist2 for some choice of dist1 and dist2. But so far I've only been talking about the probability monad. What about the quantum monad Q? In that case the notion is almost identical, except the language is inverted. In that case, a distribution on (a,b) is said to be entangled if it can't be written as tensor state1 state2. Despite all of the mystification in the popular literature, that's all it means. So why is entanglement important? One of the key ideas is this. Suppose we have a process that acts only one one part of a combined system. In other words, consider something like this:

> apply_part f x = do
>   (a,b) <- x
>   b' <- f b
>   return (a,b')

If the two subcomponents are independent (or unentangled) so that x=tensor dist1 dist2 then apply_part f x=tensor dist1 (dist2 >>= f). In other words, if we're dealing with a distribution like tensor dist1 dist2, and we're interested in studying operations on the second component, then we can completely ignore the first component during the process.

> coin = 0.75 .* return True + 0.25 .* return False
> f x = 0.125 .* return x + 0.875 .* return (not x)
> ex1 = collect (tensor coin (coin >>= f))
> ex2 = collect (apply_part f (tensor coin coin))

ex1=ex2. Now suppose we have a combined system consisting of two parts: some subsystem, and its environment, ie. something of type (Environment,Subsystem). Then as long as the environment and the subsystem are independent, we can study the subsystem as if the environment weren't there. But now suppose information can 'leak' out from the subsystem. For example consider a process like:

> a `xor` b = a/=b
> process x = do
>   (e,s) <- x
>   let e' = e `xor` s
>   return (e',s)

Because there is some flow from s to e', there is no guarantee that process x has independent subparts. As a result, we can no longer treat s as an independent subsystem, and theorems that are true about an independent s no longer hold. For example, define:

> rotate :: Float -> 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

We know that return True >>= rotate (pi/10) >>= rotate (-pi/10) gives us back the state we started with. Similarly, with

> pure x = do
>   (a,b) <- x
>   b' <- rotate (pi/10) b
>   b'' <- rotate (-pi/10) b'
>   return (a,b'')

pure acts as the identity on the second component. In particular, pure (e,True) is the same as tensor e (return True) because the terms involving false cancel out. But if we have a 'leak' as in:

> impure x = do
>   (a,b) <- x
>   b' <- rotate (pi/10) b
>   a' <- return (a `xor` b') -- information flowing out
>   b'' <- rotate (-pi/10) b'
>   return (a',b'')

and consider impure (e,True) we end up with something where the False terms don't cancel out. Compare collect $ pure $ return (False,False) with collect $ impure $ return (False,False). We're familiar with the idea that something flowing into a system might mess it up. But here we have the situation that information flowing out of a physical system can also mess it up. Any kind of entanglement between a system and its environment is called decoherence. And what I've shown is that interactions between a system and its environment can cause a physical system to 'decohere', even if the direction of interaction appears to be only outwards.

And that's one reason why quantum computing is so hard. Not only must nothing in the environment affect our qubits, neither must any information about our qubits leave the system before we're ready to have them do so.

Anyway, now I'm equipped to start writing code to simulate quantum error correction...

Update: Tried to reuse some of this code and found a bug which I've now fixed. Sorry if this pops it to the front of various feeds. It's an old post.

Monads, Vector Spaces and Quantum Mechanics pt. II

Back from wordpress.com:

I had originally intended to write some code to simulate quantum computers and implement some quantum algorithms. I'll probably eventually do that but today I just want to look at quantum mechanics in its own right as a kind of generalisation of probability. This is probably going to be the most incomprehensible post I've written in this blog. On the other hand, even though I eventually talk about the philosophy of quantum mechanics, there's some Haskell code to play with at every stage, and the codeives the same results as appear in physics papers, so maybe that will help give a handle on what I'm saying.

First get some Haskell fluff out of the way:

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

Now define certain types of vector spaces. The idea is the a W b a is a vector in a space whose basis elements are labelled by objects of type a and where the coefficients are of type p.

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

This is very similar to standard probability monads except that I've allowed the probabilities to be types other than Float. Now we need a couple of ways to operate on these vectors.

mapW allows the application of a function transforming the probabilities...

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

and fmap applies a function to the basis element labels.

> instance Functor (W b) where

>    fmap f (W a) = W $ map (\(a,p) -> (f a,p)) a

We want our vectors to support addition, multiplication, and actually form a monad. The definition of >>= is similar to that for other probability monads. Note how vector addition just concatenates our lists of probabilities. The problem with this is that if we have a vector like a+2a we'd like it to be reduced to 3a but in order to do that we need to be able to spot that the two terms a and 2a both contain multiples of the same vector, and to do that we need the fact that the labels are instances of Eq. Unfortunately we can't do this conveniently in Haskell because of the lack of restricted datatypes and so to collect similar terms we need to use a separate collect function:

> 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 . toList . fromListWith (+) . runW

Now we can specialise to the two monads that interest us:

> type P a = W Float a

> type Q a = W (Complex Float) a

P is the (hopefully familiar if you've read Eric's recent posts) probability monad. But Q allows complex probabilities. This is because quantum mechanics is a lot like probability theory with complex numbers and many of the rules of probability theory carry over.

Suppose we have a (non-quantum macroscopic) coin that we toss. It's state might be described by:

> data Coin = Heads | Tails deriving (Eq,Show,Ord)

> coin1 = 0.5 .* return Heads + 0.5 .* return Tails :: P Coin

Suppose that if Albert sees a coin that is heads up he has a 50% chance of turning it over and if he sees a coin that is tails up he has a 25% chance of turning it over. We can describe Albert like this:

> albert Heads = 0.5 .* return Heads + 0.5 .* return Tails
> albert Tails = 0.25 .* return Heads + 0.75 .* return Tails

We can now ask what happens if Albert sees a coin originally turned up heads n times in a row:

> repeat 0 f = id
> repeat n f = repeat (n-1) f . f
> (->-) :: a -> (a -> b) -> b

> g ->- f = f g

> (-><) :: Q a -> (a -> Q b) -> Q b
> g ->< f = g >>= f

> albert1 n = return Heads ->- repeat n (->< albert) ->- collect

Let me explain those new operators. ->- is just function application written from left to right. The > in the middle is intended to suggest the direction of data flow. ->< is just >>= but I've written it this way with the final < intended to suggest the way a function a -> M b 'fans out'. Anyway, apropos of nothing else, notice how Albert approaches a steady state as n gets larger.

Quantum mechanics works similarly but with the following twist. When we come to observe the state of a quantum system it undergoes the following radical change:

> observe :: Ord a => Q a -> P a

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

Ie. the quantum state becomes an ordinary probabilistic one. This is called wavefunction collapse. Before collapse, the complex weights are called 'amplitudes' rather than probabilities. The business of physicists is largely about determining what these amplitudes are. For example, the well known Schrödinger equations is a lot like a kind of probabilistic diffusion, like a random walk, except with complex probabilities instead of amplitudes. (That's why so many physicists have been hired into finance firms in recent years - stocks follow a random walk which has formal similarities to quantum physics.)

The rules of quantum mechanics are a bit like those of probability theory. In probability theory the sum of the probabilites must add to one. In addition, any process (like albert) must act in such a way that if the input sum of probabilities is one, then so is the output. This means that probabilistic process are stochastic. In quantum mechanics the sum of the squares of the magnitudes of the amplitudes must be one. Such a state is called 'normalised'. All processes must be such that normalised inputs go to normalised outputs. Such processes are called unitary ones.

There's a curious subtlety present in quantum mechanics. In classical probability theory you need to have the sum of the probabilities of your different events to sum to one. But it's no good having events like "die turns up 1", "die turns up 2", "die turns up even" at the same time. "die turns up even" includes "die turns up 2". So you always need to work with a mutually exclusive set of events. In quantum mechanics it can be pretty tricky to figure out what the mutually exclusive events are. For example, when considering the spin of an electron, there are no more mutually exclusive events beyond "spin up" and "spin down". You might think "what about spin left?". That's just a mixture of spin up and spin down - and that fact is highly non-trivial and non-obvious. But I don't want to discuss that now and it won't affect the kinds of things I'm considering below.

So here's an example of a quantum process a bit like albert above. For any angle $latex \theta$, rotate turns a boolean state into a mixture of boolean states. For $latex \theta=0$ it just leaves the state unchanged and for $latex \theta=\pi$ it inverts the state so it corresponds to the function Not. But for $latex \theta=\pi/2$ it does something really neat: it is a kind of square root of Not. Let's see it in action:

> rotate :: Float -> 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

> snot = rotate (pi/2)

> repeatM n f = repeat n (>>= f)

> snot1 n = return True ->- repeatM n snot ->- observe

We can test it by running snot1 2 to see that two applications take you to where you started but that snot1 1 gives you a 50/50 chance of finding True or False. Nothing like this is possible with classical probability theory and it can only happen because complex numbers can 'cancel each other out'. This is what is known as 'destructive interference'. In classical probability theory you only get constructive interference because probabilities are always positive real numbers. (Note that repeatM is just a monadic version of repeat - we could have used it to simplify albert1 above so there's nothing specifically quantum about it.)

Now for two more combinators:

> (=>=) :: P a -> (a -> b) -> P b
> g =>= f = fmap f g
> (=><) :: P (Q a) -> (a -> Q b) -> P (Q b)
> g =>< f = fmap (>>= f) g

The first just uses fmap to apply the function. I'm using the = sign as a convention that the function is to be applied not at the top level but one level down within the datastructure. The second is simply a monadic version of the first. The reason we need the latter is that we're going to have systems that have both kinds of uncertainty - classical probabilistic uncertainty as well as quantum uncertainty. We'll also want to use the fact that P is a monad to convert doubly uncertain events to singly uncertain ones. That's what join does:

> join :: P (P a) -> P a
> join = (>>= id)

OK, that's enough ground work. Let's investigate a physical process that can be studied in the lab: the Quantum Zeno effect, otherwise known as the fact that a watched pot never boils. First an example related to snot1:

> zeno1 n = return True ->- repeatM n (rotate (pi/fromInteger n)) ->- collect ->- observe

The idea is that we 'rotate' our system through an angle $latex \pi/n$ but we do so in n stages. The fact that we do it in n stages makes no difference, we get the same result as doing it in one go. The slight complication is this: suppose we start with a probabilistic state of type P a. If we let it evolve quantum mechanically it'll turn into something of type P (Q a). On observation we get something of type P (P a). We need join to get a single probability distribution of type P a. The join is nothing mysterious, it just combines the outcome of two successive probabilistic processes into one using the usual laws of probability.

But here's a variation on that theme. Now we carry out n stages, but after each one we observe the system causing wavefunction collapse:

> zeno2 n = return True ->- repeat n (
>     \x -> x =>= return =>< rotate (pi/fromInteger n) =>= observe ->- join
>     ) ->- collect

Notice what happens. In the former case we flipped the polarity of the input. In this case it remains closer to the original state. The higher we make n the closer it stays to its original state. (Not too high, start with small n. The code suffers from combinatorial explosion.) Here's a paper describing the actual experiment. Who needs all that messing about with sensitive equipment when you have a computer? :-)

A state of the form P (Q a) is called a mixed state. Mixed states can get a bit hairy to deal with as you have this double level of uncertainty. It can get even trickier because you can sometimes observe just part of a quantum system rather than the whole system like oberve does. This inevitably leads mixed states. von Neumann came up with the notion of a density matrix to deal with this, although a P (Q a) works fine too. I also have a hunch there is an elegant way to handle them through an object of type P (Q (Q a)) that will eliminate the whole magnitude squared thing. However, I want to look at the quantum Zeno effect in another way that ultimately allows you deal with mixed states in another way. Unfortunately I don't have time to explain this elimination today, but we can look at the general approach.

In this version I'm going to consider a quantum system that consists of the logical state in the Zeno examples, but also include the state of the observer. Now standard dogma says you can't can't form quantum states out of observers. In other words, you can't form Q Observer where Observer is the state of the observer. It says you can only form P Observer. Whatever. I'm going to represent an experimenter using a list that representing the sequence of measurements they have made. Represent the complete system by a pair of type ([Bool],Bool). The first element of the pair is the experimenter's memory and the second element is the state of the boolean variable being studied. When our experimenter makes a measurement of the boolean variable, its value is simply prepended to his or her memory:

> zeno3 n = return ([],True) ->- repeatM n (
>    \(m,s) -> do
>     s' <- rotate (pi/fromInteger n) s
>     return (s:m,s')
>    ) ->- observe =>= snd ->- collect

Note how we now delay the final observation until the end when we observe both the experimenter and the poor boolean being experimented on. We want to know the probabilities for the final boolean state so we apply snd so as to discard the state of the observer's memory. Note how we get the same result as zeno2. (Note no mixed state, just an expanded quantum state that collapses to a classical probabilistic state.)

There's an interesting philosophical implication in this. If we model the environment (in this case the experimenter is part of that environment) as part of a quantum system, we don't need all the intermediate wavefunction collapses, just the final one at the end. So are the intermediate collapses real or not? The interaction with the environment is known as decoherence and some hope that wavefunction collapse can be explained away in terms of it.

Anyway, time for you to go and do something down-to-earth like gardening. Me, I'm washing the kitchen floor...

---

I must mention an important cheat I made above. When I model the experimenter's memory as a list I'm copying the state of the measured experiment into a list. But you can't simply copy data into a quantum register. One way to see this is that unitary processes are always invertible. Copy data into a register destroys the value that was there before and hence is not invertible. So instead, imagine that we really have an array that starts out zeroed out and that each time something is added to the list, the new result is xored into the next slot in the array. The list is just non-unitary looking, but convenient, shorthand for this unitary process.

The allure of LaTeX

I'm moving! Get over here. I don't know if this will be a permanent move but I'll try it for a bit. Old posts haven't been imported too well so I may just keep the old posts here and delete the ones I imported there.

Monads in C pt. III

OK, the last C monad example:

typedef struct {
  char *string;
  int value;
} WriterInt;

WriterInt returnWriter(int i)
{
  WriterInt r;
  r.string = "";
  r.value = i;
  return r;
}

WriterInt bind(WriterInt (*f)(int),WriterInt x)
{
  WriterInt y = f(x.value);
  WriterInt z;
  z.value = y.value;
  int len = strlen(x.string)+strlen(y.string);
  z.string = malloc(len+1);
  strcpy(z.string,x.string);
  strcat(z.string,y.string);
  return z;
}

WriterInt print(int i)
{
  WriterInt x;
  x.string = malloc(32);
  x.value = sprintf(x.string,"%d\n",i);
  return x;
}

WriterInt printplus_bad(int i)
{
  WriterInt x = print(i);
  return print(x.value); /* cheating! */
}

WriterInt printplus(int i)
{
  WriterInt x = print(i);
  return bind(print,x);
}

This time, instead of printing, we build up a string as a side effect. printplus() is implemented exactly as before, without knowledge of how to handle side effects, and yet it correctly handles strings returned as a side effect and concatenates them together. The magic is bind() which allows a function expecting an int input to be applied to a WriterInt.

I'm hoping that programmers of languages such as C can now begin to see the pattern shared by these examples.

Anyway, this is all just stalling for time while I try to get my quantum computation code working.

Monads in C, pt. II

This is a quick followup to a post on reddit. This time a crude approximation to the Maybe monad implemented in C:

include <stdio.h>

typedef struct {
   int something;
   int just;
} MaybeInt;

MaybeInt returnMaybe(int i)
{
   MaybeInt r;
   r.something = 1;
   r.just = i;
   return r;
}

MaybeInt nothing()
{
   MaybeInt r;
   r.something = 0;
   return r;
}

MaybeInt bind(MaybeInt (*f)(int),MaybeInt x)
{
   if (x.something)
   {
       return f(x.just);
   } else
   {
       return nothing();
   }
}

MaybeInt print(int i)
{
   int written = printf("%d\n",i);
   if (written>=0)
   {
       return returnMaybe(written);
   } else
   {
       return nothing();
   }
}

MaybeInt printplus_bad(int i)
{
   MaybeInt x = print(i);
   return print(x.just); /* cheating! */
}

MaybeInt printplus(int i)
{
   MaybeInt x = print(i);
   return bind(print,x);
}

Again the idea is that printplus() is a version of printplus_bad() that uses just the two function bind/return interface to the MaybeInt type to achieve its effect. This time, instead of simply tainting IO with a particular type, the Maybe monad is able to deal gracefully with failure. If, for some obscure reason, printf() fails, returning an integer less than zero, it returns an object representing this fact. If, as in printplus(), you have two calls to print() in a row, bind() handles all the plumbing for you automatically. The really important thing is this: the implementation of printplus() is identical to my previous example, and yet the semantics is quite different because printplus() is able to bail out early. This bailing out is completely hidden inside bind().

I hope that this gives some hint of what monads can do to even hardcore non-functional programmers. If not, I'll probably write another example soon.

(Remember of course that this isn't meant to be practical code. It was a response to someone who wanted to at least some some C code for monads to get an idea of what they're about.)