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

Monads for vector spaces, probability and quantum mechanics pt. I

I've been enjoying Eric Kidd's articles on probability theory with Haskell. So I thought I'd follow them up with two things: (1) finding the underlying algebraic structure and (2) showing how it's general enough to do more than just probability.

Firstly, Eric found a really neat factoring of the probability monad as a 'product' of two monads: the Perhaps monad and the List monad. This factoring can be interpreted algebraically.

A monoid is defined to be a set m with a binary operator (typically written as abuttal, ie. the product of a and b is ab) that is associative, and with an element, 1, that is an identity for the binary operator. The List monad gives rise to the monoid freely generated by a set, in other words it defines the smallest monoid containing the set and with no extra relationships between the elements that don't come from the definition of a monoid. The binary operator for List is called ++ and we embed the original set in the monoid using the function \x -> [x], otherwise known as return. [] is the identity. It's not hard to see that we have asociativity as (a ++ b) ++ c == a ++ (b ++ c). If we use Set instead of List then we get the free commutative monoid instead, ie. one where a+b=b+a.

(Eh? Where's my paragraph on m-sets. Don't tell me I have to format all those HTML subscripts again! I hate it when computers do that. I wrote a paragraph and it just vanished. Oh well, maybe version 2.0 will be better. But without subscripts this time...)

If m is a monoid then an m-set is a set s with an action of m on it. An action of a monoid on a set is a scheme that converts each element of the monoid, say g, into a map g^*:S→S so that 1^* is the identity map and g^*h^*=(gh)^*. The type Writer m s corresponds to pairs in s×m but this also doubles as the free m-set. We simply define f^*(x,g)=(x,fg). It's free because we never get f^*(x,f')=g^*(x',g') unless x=x', so we never get any surprising equalities that aren't inherent in the definition of an m-set. Perhaps, which is actually a pseudonym for Writer Float, defines an R-set, where R is the monoid of reals under multiplication.

In an earlier post I showed how you could layer up algebraic structures. There I showed that you could combine a commutative monoid structure with a monoid structure to get a semiring. So what happens if we layer up a commutative monoid and an m-set? We'll get something that has a commutative binary operator but that can also be acted on my elements of m. Specialise to the case when m is a field like the real or complex numbers. Does this sound like a familiar structure? I hope so, it's a description of a vector space. What Eric Kidd has shown is that we can build a vector space monad by applying an m-set monad transformer (with m a field) to a commutative monoid monad. (This isn't a completely trivial fact and it depends on the fact that the definition of PerhapsT handles distributivity correctly.) I think that's pretty neat. But I'm getting slightly ahead of myself here as I haven't shown what Eric's stuff has to do with vector spaces.

A probability distribution can be thought of as a vector with outcomes forming a basis. Any distribution attaches a probability 'weight' to each possible outcome in the same way that a vector can be written as a weighted sum of basis elements. Each time we have a probabilistic transition, we're effectively multiplying our distribution vector by a stochastic matrix. Eric (and others)'s monad allows you to write these matrix multiplications in a very natural way.

A vector space forms a monad in a straightforward way. If V(B) is the vector space generated by basis B, and V(C) is another vector space, then any function B→V(C) can be lifted to a linear map V(B)→V(C). This lift is clearly of type (B→V(C))→(V(B)→V(C)). With a little flip we get V(B)→(B→V(C))→V(C). This is the main step in showing V to be a monad. And what the probability monad allows us to do is write our stochastic matrices in terms of what happens to the individual basis elements (ie. outcomes) instead of having to write out the entire matrix.

Anyway, this is all just waffle and it really needs some code to make it more concrete. I have an ulterior motive here. It's not just probability theory that is described by vector spaces. So is quantum mechanics. And I plan to work my way up to defining quantum computers and implementing a bunch of well known quantum algorithms in Haskell. In the next installment I expect to get at least as far as writing the code to play with the square root of NOT operator.

PS This idea of layering up algebraic structures is one I haven't found in the textbooks yet. I'm guessing it's in Mac Lane, but I haven't yet justified the cost of that book to myself. But I don't actually have a clue what a category theorist would call a monad transformer and don't recall reading about them in any texts that weren't specifically about computer science. Maybe someone else can fill me in. I do know that it's absolutely standard to get free algebraic structures from monads.

Errata: As notfancy points out, I should really be talking about multisets rather than sets, otherwise we have a+a=a. As programmers often use lists to represent both sets and multisets it can sometimes get confusing.

The Essence of Quantum Computing

A popular view of quantum computation is that its power (at least the power of certain well known quantum algorithms) derives from a form of parallel computation coming from the fact that the state of a quantum computer can be a superposition of many classical states. I want to argue that this is in fact a bit of a red herring because there are other models of computation that allow formal superpositions and yet clearly don't provide a form of parallel computation. Instead, the essence of quantum computing is derived from a kind of non-convexity property and this is the common feature of all of the well known quantum algorithms that isn't shared by classical algorithms.

Before continuing there's something important I should say. Let me quote Bennett et al. in Strengths and Weaknesses of Quantum Computing

"The power of quantum computers comes from their ability to follow a coherent superposition of computation paths."

In other words, what I am saying is in direct contradiction with what is said by some of the founding fathers of quantum computing. Actually, I think that's a good thing, it means that whether I'm right or wrong, I must be saying something non-trivial.

Consider a probabilistic computer. We can represent the probability of its being in any state by an element of a vector space whose basis is the set of possible states and where the coefficient of each basis element is the probability of the computer being in that state. Given any pair of vectors that represent possible ditributions of states. v₁ and v₂, then av₁+bv₂, where a,b>=0 and a+b = 1 is also a possible distribution.In other words we can form superpositions exactly analogously to the way we can form superpositions of quantum computer states. The mathematics that describes the notion of a quantum computer performing computations in parallel by superposing N states is formally almost identical to the mathematics describing a probabilistic computer as a superposition of states.

So if we can superpose N states in both classical probabilistic computers and quantum computers, why do neither seem to offer an N times speedup? In the case of quantum computers one popular explanation is that then though we run N states in parallel, we can't perform I/O with those N states in parallel. When we observe a quantum computer we cause a wavefunction collapse (it doesn't matter whether you are a believer in such collapse, that is how we model quantum computers) and instead we make a single observation from a probability distribution. Mathematically this is identical to what happens in a probabilistic machine - a superposition is nondeterministically mapped into one of many definite states.

So whatever the difference is between probabilistic computers and quantum computers, it doesn't lie in the fact that you can form superpositions. Mathematically speaking, probabilistic computers and quantum computers are almost identical in this respect.

But we know that quantum algorithms can do things that no probabilistic algorithm can, so there must be some difference. There is of course. Here's an approach to seeing what it is.

Consider a probabilistic branching process:

    __ A
   /  
  2
 / \__ B
1   __ A
 \ / 
  3
   \__ C

Starting at state 1 we can either go to state 2 or state 3, and from there we end up at either A, B or C. Assume that the choice of branch is nondeterministic and that we can assign probabilities to each branch. Suppose also that the branch probabilities are all non-zero. The important thing to note in that diagram is that there are two ways to reach it. The probability of reaching A is the sum of the probabilities for each of the different ways to reach A. The key point I want to make here is that A is more likely because there are two ways to reach it. If I were to replace one or other of the A's with a D then the likelihood of finding our process in state A at the end is reduced.

But in quantum mechanics we replace probabilities with complex 'amplitudes'. When we make our observation at the end, the probability of any outcome is given by the modulus of the amplitude squared. The amplitudes still add, just as in probability theory. But this has a very important consequence. If the amplitudes for the two different ways to reach A lie in opposite directions in the complex plane, then the probability of ending up in state A after observation is reduced by having two ways of reaching A. This is completely different to the probabilistic case. We can paradoxically make something less likely by adding a new way to achieve it. This is where the essence of quantum computing lies.

Pick up just about any paper describing a quantum algorithm and you'll see that the important step is a 'conspiracy' to make the amplitudes of 'undesirable' states sum to zero (or at least something small). Probably the simplest example is counterfactual computation where it is arranged that the computation is executed along two different branches with equal and opposite amplitudes. The neatest example is probably Shor's cycle discovery algorithm that forms the basis of his factoring algorithm. By clever use of the Fourier transform he arranges that branches leading to non-cycles cancel each other out leaving only the periods of cycles with high amplitudes.

However, I think it's important to stress that given any individual quantum computer Bennett et al. and I probably don't disagree on what it actually does. This is an issue of semantics. But semantics does make a difference. The notion of a quantum computer as a super-parallel machine is one being promulgated in the popular science press and I think this causes much confusion with people expecting all kinds of algorithms to run in O(1) time.

So why do people talk about parallel computation if that's not what quantum computer do? I think it's because people aren't comfortable with quantum computers yet. This means they don't want to think of quantum computers as something fundamental and irreducible and instead want to reduce them to something classical. The easiest classical model is the parallel one and so they 'reduce' it to that.

PS This is essentially a specialisation of an old post to the context of quantum computing.

Modular arithmetic with regular expressions

Problem: Find a regular expression, compatible with GNU grep, that recognises strings of 1's and 0's that form the binary expansion of multiples of 7. (Consider the empty string to be zero.)

Solution:
Here's one (though not necessarily a good one):

^\(0\|1\(1\(01*00\)*01*010\|\(0\|1\(01*00
\)*01*011\)\(10\(01*00\)*01*011\|\(0\|11\
)1\)*\(10\(01*00\)*01*010\|\(0\|11\)0\)\)
*\(1\(01*00\)*1\|\(0\|1\(01*00\)*01*011\)
\(10\(01*00\)*01*011\|\(0\|11\)1\)*10\(01
*00\)*1\)\)*$

The standard approach to solving this problem is probably to construct a finite state automaton to compute the remainder modulo 7 of a binary string and then use the proof of Kleene's theorem to convert the automaton to a regular expression. But as I have some regular expression code left over from a previous post I thought I'd try a purely algebraic approach. Although it's ultimately equivalent to deriving a finite state automaton I thought it might be more fun to use intuitions about linear algebra to get to the solution. And anyway, drawing diagrams with HTML is hard...

Consider first recognising multiples of 3. Assume we have three regular expressions, A, B and C. A corresponds to multiples of 3, B corresponds to numbers that are 1 modulo 3 and C matches numbers that are 2 modulo 3. Then we know

A = A0+B1+1

B = C0+A1

C = B0+C1

First a note about notation: I'm using bold 0 and 1 to represent binary digits but using regular 0 to represent the regular expression matching nothing and regular 1 to represent the regular expression matching the empty string. a+b means match a or b. So we can interpret the equation for A, say, as saying that a multiple of 3 is the empty string, or a multiple of 3 followed by a 0, or a 1 after one more than a multiple of 3. Similar interpretations apply for the other two equations. We can rewrite this as

x = xM+v

where

M=
0 1 0
1 0 0
0 0 1

and

x = [A B C], v = [1 0 0]

We can now use our ordinary intuitions about solving linear systems. We want to solve for A so we need to eliminate B and C. Start by eliminating C. We have

C = B0+C1

If we were doing ordinary linear algebra we could say

C(1-1) = B0

C = B0(1-1)^-1

but for regular expressions we can neither subtract nor divide. However, arguing informally now, (1-a)^-1=1+a+a²+a³+…=a^* so now we can write

C = B01^*

Anyway, those are all the principles we need. Just eliminate each variable in turn by solving for it in terms of the other remaining variables and substituting the solution back into the remaining equations. You should be able to finish the job and solve for A. But you'll notice the work gets tedious for larger integers. That's where the code comes in.

It starts as before but we can delete a whole lot of material:

> import List
> data RE a = Symbol a | Star (RE a) | Sum [RE a] | Prod [RE a]
>   deriving (Eq,Show,Ord)

> a <*> b = Prod [a,b]
> a <+> b = Sum [a,b]

> zero = Sum []
> one = Prod []

My simplification rules are still really poor. (Buggy, in fact, but not buggy in a way that leads to incorrect results.) That means we don't expect the final regular expressions to be in any way efficient. But I'm not trying to get the best solution, just a solution.

> simplify a = let b = simplify' a in if a==b then b else simplify b where
>   simplify' (Prod [a]) = simplify a
>   simplify' (Prod a) | zero `elem` a = zero
>   simplify' (Prod (a:b)) = case simplify a of
>       Prod x -> Prod $ x ++ map simplify b
>       a -> Prod $ filter (/= one) $ map simplify (a:b)

>   simplify' (Sum [a]) = a
>   simplify' (Sum (a:b)) = case simplify a of
>       Sum x -> Sum $ x ++ map simplify b
>       a -> Sum $ nub $ sort $ filter (/= zero) $ map simplify (a:b)

>   simplify' (Star (Sum [])) = one
>   simplify' (Star (Prod [])) = one
>   simplify' (Star a) = Star (simplify a)
>   simplify' a = a

I'm representing the vector v and matrix M above as simple associative arrays whose elements are found using lookup. Elements not represented are assumed to be 0. fetch handles the sparse lookup:

> type Matrix s a = [((s,s),RE a)]
> type Vector s a = [(s,RE a)]

> fetch i x = maybe zero id $ lookup i x

solve implements exactly the solution method described above. The use of Star corresponds to my abuse of the binomial series for (1-a)^-1 above. However rather than manipulate equations I'm always working in terms of matrices and vectors over the regular expressions:

> solve range i (m,v) = (
>  [(j,simplify $ fetch (i,j) m <*> Star (fetch (i,i) m)) |
>        j <- range, j/=i],
>        simplify $ fetch i v <*> Star (fetch (i,i) m)
>     )

substituteFor substitutes a solution to an equation given by solve into the remainder of the equations.

> substituteFor range k (v',c') m v = (
>  [((i,j),(fetch j v' <*> fetch (i,k) m) <+> fetch (i,j) m) |
>         i <- range, i/=k, j <- range, j/=k],
>  [(i,(c' <*> fetch (i,k) m) <+> fetch i v) |
>         i <- range, i /= k]
>     )

As I describe above, a variable is eliminated by solving for it and then substituting:

> elim range s (m,v) = let (v',c') = solve range s (m,v)
>                in substituteFor range s (v',c') m v

Now some code to output our regular expressions in a form compatible with grep:

> bracket n (s,n') = if n'>=n then s else "\\(" ++ s ++ "\\)"
> pretty (Symbol a) = ([a],9)
> pretty (Star a) = (bracket 8 (pretty a) ++ "*",8)
> pretty (Prod a) = (concat $ map (bracket 7) $ map pretty a,7)
> pretty (Sum a) = (foldr1 (\a b -> a ++ "\\|" ++ b) $
>                   map (bracket 6) $ map pretty a,6)

Now construct M and v for the equations describing the regular expressions for the classes of binary strings modulo n:

> m n = [(((2*i) `mod` n,i),Symbol '0')|i <- [0..n-1]] ++
>       [(((2*i+1) `mod` n,i),Symbol '1')|i <- [0..n-1]]
> v n = [(0,one)]

To solve a complete set of equations we reduce it one equation at a time by eliminating variables:

> reduce [a] (m,v) = solve [a] a (m,v)
> reduce (a:as) (m,v) = reduce as $ elim (a:as) a (m,v)
> test n = let (_,r) = reduce [n-1,n-2..0] (m n,v n) in simplify r

Print regular expression to detect multiples of 7:

> main = putStrLn $ "^" ++ fst (pretty (test 7)) ++ "$"

Try integers besides 7. But note that the expressions get more and more complex. And of course you can use the equation solving approach for problems besides recognising multiples of integers.

PS I just noticed some of the indentation went wrong on translation to HTML. But the code appears to still compile fine with GHC.

Exceptions, Disjunctions and Continuations

At the end of section 4 of Call-by-value is dual to call-by-name, Wadler tells a story about the Devil and a billion dollars. After a long time thinking about it, the meaning of it suddenly hit me like an epiphany a few weeks ago, and I'm finally in a position to give an explanation of what it means. (This is also an exercise in TMR6.) And just for a bit of variety I'm going to throw in some C++ and Scheme, starting with some simple C++ code:

#include <iostream>
#include <string>

using namespace std;

template<class A,class B>
class Either
{
public:
    virtual B value() const throw (A) = 0;
};

template<class A,class B>
class Left : public Either<A,B>
{
    A a;
public:
    Left(A a0) : a(a0) { }
    B value() const throw (A) { throw a; }
};

template<class A,class B>
class Right : public Either<A,B>
{
    B b;
public:
    Right(B b0) : b(b0) { }
    B value() const throw (A) { return b; }
};

template<class A,class B,class C>
C either(const Either<A,B> &e,C (*f)(A),C (*g)(B))
{
    try {
        return g(e.value());
    } catch (A &a) {
        return f(a);
    }
}

string f(int a)
{
    return "int";
}

string g(float b)
{
    return "float";
}

int main()
{
    std::cout <<
        either<int,float,string>(Left<int,float>(7),f,g)
        << std::endl;
    std::cout <<
        either<int,float,string>(Right<int,float>(7.0f),f,g)
        << std::endl;
}

This provides a (somewhat unreliable, I think) implementation of a discriminated union. An object of type Either<A,B> is either an object of type A or an object of type B. But unlike the ordinary union, it remembers the type of its contents. I've also provided an either function. Given a function of type A→C and a function of type B→C, then if we're handed an object of type Either<A,B> we know that we can apply one or other function to it to get an object of type C. The functon either applies the appropriate one for you. What's unusual is the way I've implemented it. I'm representing an object of type Either<A,B> as a function of no arguments (in effect a lazy value) returning something of type B (*)() throw (A). There are two paths by which an object can get out of such a function: either it is returned or it is thrown out. Having two paths is what I've used to implement of a template that represents either one thing or another. This simple example should make it clear that there is some sort of a connection between exceptions and logical or (ie. disjunction).

In functional languages exceptions tend to be handled a little differently. Whereas C++ has a throw keyword, there is no equivalent in Scheme, say. Instead
the call-with-current-continuation function acts like catch. The argument to call-with-current-continuation is a function (lambda (throw) ...) and the argument to this function (itself a function), here called throw, does the throwing. For example consider the scheme expression:

(call-with-current-continuation
    (lambda (throw) (+ 1 (throw 7))))

The result is 7, not 8. The 7 is thrown 'past' the (+ 1 ...) and is caught by the surrounding call-with-current-continuation. The important difference from C++ is that if any subexpression is to throw some value, it must be 'handed' a suitable throw function. Note also that Scheme is dynamically typed. In a static functional language throwing exceptions is a tiny bit more subtle. Not only must a subexpression be handed a suitable function for throwing, it must be handed a function of a suitable type. In particular, if the data thrown is to be of type a, then the subexpression must be handed a suitable 'throw' function that can take an object of type a as argument. Think of such a function as an a-consumer. So if we want to implement the C++ idea above in a statically typed functional language that supports exceptions we need to replace the function of type B (*)() throw (A) with a function that consumes an A-consumer and returns an object of type B. Such a function can either throw a result of type A or return a result of type B.

Anyway, here's a complete Scheme example (BTW it's my first bit of Scheme code longer than one line):

(define (right b) (lambda (throw) b))
(define (left a) (lambda (throw) (throw a)))

(define (compose f g) (lambda (x) (f (g x))))
(define (id x) x)

(define (either f g x)
   (call-with-current-continuation
       (lambda (throw)
           (g (x (compose throw f))))))

(define (f x) (cons "left" x))
(define (g x) (cons "right" x))

(either f g (left 1))
(either f g (right 2))

OK, time for some logic.

As everyone familiar with propositional calculus knows, a ∨ b = ~a → b. Or in English "a or b" is the same as "if a isn't true then b is".

Less well known is the Curry-Howard isomorphism. In a nutshell it goes like this: if you know a, and you know a → b, then you can deduce b (using modus ponens). Similarly, in a strictly typed programming language, if you have an object of type a, and a function of type a → b, then you can use them to make an object of type b (just apply the function to a). So given a bunch of propositions in 1-1 correspondence with a bunch of objects, you should be able to see that the set of propositions you can deduce by modus ponens is going to be in 1-1 correspondence with the objects you can make by function application.

But there's more to logic than modus ponens. In classical logic you also have double negation elimination. In other words, if you know ~~a then you can deduce a. So double negation elimination is like having a polymorphic function ~~a → a for any a. What does "~" mean in terms of functions? In logic, define ⊥ to be the unprovable false proposition. Then ~a is simply a→⊥. In programming define ⊥ to be the type with no instances and we can define ~a to be the type of functions a→⊥. A function of type ~a is a function that consumes an instance of type a and gives you back nothing in return. Clearly if you can make an object of type a then no function of type ~a can exist or else you could apply it to a and get an object f type ⊥, which is impossible by definition.

It's not hard to see that there can be no well behaved function that for any type a can map an ~~a to an a. An ~~a is a thing that consumes a-consumers. In a sense it's destructive. There are no raw ingredients from which you could possibly make an a for all a. So there is no function of type ~~a → a. But this type can serve very well as the type of the catch operator described above. Rewriting the scheme example above with Haskell-like notation we expect something like

catch $ \throw -> ... (throw a) ...

Because throw never returns we can assign it the type a -> ⊥ without problems and hence catch can be assigned type ~~a. This now means that the Curry-Howard isomorphism can be extended to classical logic if we throw in catch along with the usual functions you expect in a pure typed functional language.

So If we use catch to throw an exception out of a function, then the function must be handed a ~a as argument. So in order to implement Either, we must replace the C++ type B (*)() throw(A) with ~a -> b. In other words, if we use the symbol ∨ to mean either, then exactly in classical logic, we have a ∨ b = ~a → b.

You can think of the argument to catch as being a bit like a bungee cord with one end attached to the catch. If a function is handed the other end of such a bungee cord it can either use it to catapult out a value of type a, or the function can ignore it and return somehing of type b. But there's a neat twist on this. If b = ~a then the function, on being handed a bungee cord, could simply hand it back as a return value. In other words, whatever type a is, we can always construct an object of type a ∨ ~a from the identity function. This corresponds to the fact that a ∨ ~a is a tautology of classical logic.

And now you can reread Wadler's story about the Devil and see that my story about a bungee cord is really the same one.

Exercise: Implement Left, Right and either using the method above in the Haskell Cont monad.

PS Reasoning about exceptions in this way is new to me and so I don't feel I have the above on a firm footing. So expect some errors. In particular I haven't thought about what happens if you try to nest disjunctive structures this way. You can be sure that the C++ implementation will fail because the exception may be caught by the wrong try.

Update: fixed some typos. But also here's an exercise:

Exercise: a ∧ b = ~(~a ∨ ~b). What does that mean in the above context? (I haven't got my head around this yet.)

Comonads and reading from the future

I'm too busy reading Bernie Pope and Russell O'Connor's articles in the latest Monad Reader to think about much else right now, so I'm only going to make a brief post. I want to mention something pointed out by Nick Frisby over here. The loeb function is in fact closely related to cfix, the dual of mfix for monads, defined by Dave Menendez. (Nick Frisby's statement isn't quite correct but it's interest value outweighs its truthiness.) In fact, I can borrow Russell's language to say that whereas mfix sends data back in time, cfix and loeb read data from the future!

Think of it like this: when you use cobind with a comonad you provide some function that distills down a comonadic data structure into a single value and then cobind applies this distillation repeatedly at every point in the comonadic structure to give you another comonadic structure. My canonical example is the cellular automaton where you give a rule that defines how to compute each individual cell as a function of the whole grid, and cobind applies this throughout an entire input grid to give you back an entire output grid. The functions cfix and loeb allow you to give rules for each cell where you're allowed to make reference to the entire output grid. If you think of the output grid as something in the future that you haven't yet computed, then cfix and loeb allow you to read from the future.

Here's an example of cfix it in use (where I make [] a comonad by interpreting it as a kind of one-sided zipper):

> instance Show (x -> a)
> instance Eq (x -> a)

> instance (Num a,Eq a) => Num (x -> a) where
>     fromInteger = const . fromInteger
>     f + g = \x -> f x + g x
>     f * g = \x -> f x * g x
>     negate = (negate .)
>     abs = (abs .)
>     signum = (signum .)

> class Comonad w where
>    coreturn :: w a -> a
>    cobind :: (w a -> b) -> w a -> w b

> instance Comonad [] where
>     coreturn (x:xs) = x
>     cobind f [] = []
>     cobind f (x:xs) = f (x:xs) : cobind f xs

> cfix :: Comonad d => d (d a -> a) -> a
> cfix d = coreturn d (cobind cfix d)

> ouroboros = [2*head.tail,1+head.tail,17]
> test = cfix ouroboros

In the expression test we can use head and tail to refer to the sublist to the right of each element from within the list itself. For example, in the subexpression 2*head.tail the head.tail is being applied to cfix [1+head.tail,17].

Note how in this example (1) loeb gives you back an entire list, not just one element and (2) cfix respects the comonadic structure in the sense that you can only refer to elements to the right of the current one whereas loeb lets you refer to the entire output list. I think this helps to make clear the differences and similarities between cfix and loeb.

That wasn't too brief. But I did cheat by copying and pasting code from an old post. Anyway, back to tinkering with my implementation of dropWhile...

The Monads Hidden Behind Every Zipper

Uustalu points out that behind every zipper lies a comonad. I used this in my cellular automaton example. What he doesn't mention is that there is also a monad lurking behind every zipper and that it arises in a natural way. The catch is that it's a monad in the wrong category making it tricky to express in Haskell.

Firstly, what does "wrong category" mean? In the category of vector spaces, every object is equipped with a commutative binary + operator. We can make no such assumption in Haskell. What's more, Haskell doesn't easily allow us to restrict type constructors to instances of Num say. So we can't easily implement monads in the category of vector spaces as instances of Haskell's Monad. Nonetheless, we can always write our own bind and return constrained by Num, and that's what I'll do here. I believe there are workarounds that don't involve GADT abuse but really I just want to illustrate an idea, not write a complete working application.

As I've mentioned on #haskell, convolution is comonadic. A convolution (more specifically, a discrete convolution) is a function that maps a (doubly infinite) sequence of real numbers to another such sequence of real numbers in such a way that each number in the new sequence is a weighted sum of numbers in the old sequence and the weights depend only on the relative positions of the numbers. For example: replacing each number in the sequence with the average of it and its neighbours is a convolution.

Here's a simple example. Consider the sequence ...0,0,1,2,3,0,0,... We'll replace each number with the sum of its neighbours and twice itself. So this maps to ...,0,0+0+1,0+2*1+2,1+2*2+3,2+2*3+0,3+2*0+0,0,... The structure of this computation is very similar to the structure of cellular automaton evaluation and we can reuse that code here:

> data Zipper a = Zipper [a] [a] deriving (Eq,Show)

> left (Zipper (a:as) bs) = Zipper as (a:bs)
> right (Zipper as (b:bs)) = Zipper (b:as) bs

> class Comonad w where
>    coreturn :: w a -> a
>    cobind :: (w a -> b) -> w a -> w b

> iterate1 f x = tail $ iterate f x

> instance Comonad Zipper where
>    cobind f a = fmap f $ Zipper (iterate1 left a) (iterate right a)
>    coreturn (Zipper _ (b:_)) = b

> a = Zipper (repeat 0) ([0,1,2,3]++repeat 0)
> f (Zipper (a:_) (b:c:_)) = a+2*b+c

> test = let Zipper u v = cobind f a in take 5 v

When we run test we get [1,4,8,8,3].

Now think about the structure of this computation. We write a function f that gathers up the value at a cell and its neighbours and computes the weighted sum. cobind then applies this function all the way through our sequence to get our result. The important thing is that f pulls in the data it needs, distilling it down to one value, and that's why it has type Zipper a -> a.

There's another approach to convolution. We could push or scatter data. At each point in the sequence we could take the value and push one copy left, one copy right, and leave double the original value in place. Each cell will end up with three values left in it and all we have to do is sum the values deposited in a cell. What we want therefore is to take a value in a cell, of type a, and generate an entire zipper full of values pushed from it, ie. a Zipper a. We'd then like bind for the monad to sum over the returned zippers, staggering them right and left as appropriate. Here's a table illustrating what I mean:

0 -> 0 0 0 0 0 0
1 -> 0 1 2 1 0 0
2 -> 0 0 2 4 2 0
3 -> 0 0 0 3 6 3
0 -> 0 0 0 0 0 0
----------------
     0 1 4 8 8 3

Unfortunately, as I pointed out above, bind for a monad can't perform summations, so I have to define bind' and return' instead.

> plus (Zipper a b) (Zipper a' b') = Zipper (plus' a a') (plus' b b') where
>        plus' (a:as) (b:bs) = (a+b) : plus' as bs
>        plus' a [] = a
>        plus' [] a = a

Note that I'm going to assume that elements past the end of a list (so to speak) actually represent zero and define new versions of left and right that produce the required zeros on demand:

> left' (Zipper (a:as) bs) = Zipper as (a:bs)
> left' (Zipper [] bs) = Zipper [] (0:bs)
> right' (Zipper as (b:bs)) = Zipper (b:as) bs
> right' (Zipper as []) = Zipper (0:as) []

> tail' [] = []
> tail' a = tail a
> stagger f [] = []
> stagger f (x:xs) = x : map f (stagger f xs)
> stagger1 f x = tail' (stagger f x)

> instance Functor Zipper where
>    fmap f (Zipper a b) = Zipper (map f a) (map f b)

Here's our fake monad:

> return' a = Zipper [] [a]
> bind' f x = let Zipper a b = fmap f x
>            in foldl1 plus (stagger left' b ++ stagger1 right' a)

> a' = Zipper [] [0,1,2,3]
> f' x = Zipper [x] [2*x,x]

> test' = let Zipper u v = bind' f' a' in take 5 v

Run test' and we get the same result as before by a quite different method.

You should be able to see that in some sense comonads describe lazy computations and monads represent strict ones. In strict languages you evaluate your arguments first, and then push them into a function. In lazy languages you call the function first, and it pulls in the value of its arguments as needed. In fact, in Haskell we have a strictness monad and in ocaml you can implement a lazy comonad. By and large, it's easier to write pull algorithms in Haskell. Push algorithms are much easier to implement in a language with side effects. In fact, one of the things that Haskell beginners need to learn is how to rewrite push algorithms as pull algorithms and eliminate those side effects.

What I think is interesting here is that this dichotomy between push and pull, scatter and gather, is pervasive all the way through computing, and it's quite pretty that this deep notion it can be captured nicely by a(n almost) monad and a comonad. Just to give an illustration of how pervasive it is consider 3D rendering. We typically want to take a list of primitives and turn it into an array of pixels. If we rasterise we push the primitives, asking for each primitive, what pixels does it cover. If we ray-trace we pull on the pixels asking, for each one, what primitives cover it. Many optimisations to both of these approaches consist of making judicious decisions about which part of the pipeline should be switched from push to pull or vice versa.

BTW I have an algorithm for automatically converting some types of pull algorithm into push algorithms and vice versa. It's another one of those papers that I keep meaning to finish and submit for publication one day...

Fox's Ubiquitous Free Derivative pt. 2

As promised earlier, here are some Haskell definitions of derivatives that can be applied to regular expressions. They all give operators that might occasionally be of practical use. But what's more interesting to me is that we now have enough of an intuition pump going between types and regular expression that we can siphon off ideas from papers on types and use them with regexps.

First define a regular expression type:

> import List
> data RE a = Symbol a | Star (RE a) | Sum [RE a] | Prod [RE a]
>   deriving (Eq,Show,Ord)

I use lists instead of binary operators to represent products and sums because
addition and multiplication of regular expressions is associative.

> a <*> b = Prod [a,b]
> a <+> b = Sum [a,b]

> zero = Sum []
> one = Prod []

Here's how to make a simple regular expression from a string

> re x = Prod $ map Symbol x

Exercise: Write a proper parser for these things so we can use conventional regexp notation!

Here's a test to see whether a regular expression matches the empty string:

> acceptsEmpty (Symbol _) = False
> acceptsEmpty (Star _) = True
> acceptsEmpty (Sum l) =  or $ map acceptsEmpty l
> acceptsEmpty (Prod l) = and $ map acceptsEmpty l

I've decided to name these derivatives after people. So first I define the McBride derivative which is a generalisation of the dissection operator introduced in Clowns to the left.... The key idea is that if d = mcbride c j a then d (x <*> y) = (d x <*> j y) <+> (c x <*> d y). This needs to respect the associativity of <*> and one way to achieve this is by choosing c and j that are multiplicative homomorphisms so that, for example, c(x <*> y) = c x <*> c y.

> mcbride _ _ a (Symbol a') = if a==a' then one else zero
> mcbride c j a (Sum p) = Sum $ map (mcbride c j a) p
> mcbride _ _ _ (Prod []) = zero
> mcbride c j a (Prod x) = Sum $ map f (splits x) where
>     f (x,y:ys) = Prod $ map c x ++ [mcbride c j a y] ++ map j ys
> mcbride c j a (Star p) = c (Star p) <*> mcbride c j a p <*> j (Star p)

And now I can define all the Leibniz, Fox and Brzozowski derivatives as special cases of the mcbride derivative. The latter each come in left- and right-handed pairs.

> delta p = if acceptsEmpty p then one else zero

> leibniz,fox_l,fox_r,brzozowski_l :: (Ord a,Eq a) => a -> RE a -> RE a
> leibniz a p = simplify $ mcbride id id a p
> fox_l a p = simplify $ mcbride id (const (one)) a p
> fox_r a p = simplify $ mcbride (const (one)) id a p
> brzozowski_l a p = simplify $ mcbride delta id a p
> brzozowski_r a p = simplify $ mcbride id delta a p

brzozowski_l is the usual Brzozowski derivative. If r matches some language, brzozowski_l a t matches the sublanguage of strings beginning with a and with that a removed. We can define regular expression matching in the usual way in terms of it:

> matches re s = acceptsEmpty (matches' re s) where
>   matches' re [] = re
>   matches' re (a:as) = matches' (brzozowski_l a re) as

Now we can start matching.

ex1 would normally be written as "abc|b*c".

> ex1 = (re "abc") <+> (Star (re "b") <*> Symbol 'c')
> test1 = map (matches ex1) ["abc","bbc","bbbbc","abbbbc"]

brzozowski_r chews off a symbol from the right:

> test2 = map (matches (brzozowski_r 'c' ex1)) ["abc","bbc","ab","bbbbb"]

leibniz a allows you to leave out one a from anywhere in the regexp.

ex3 matches "abracadabra" as it might be shouted by a wizard who might lengthen the second last 'a' to make sure the spirits can hear. ex4 matches such strings with precisely one missing 'r'.

ex3 would normally be written as "abracada*bra".

> ex3 = re "abracad" <*> Star (re "a") <*> re "bra"
> test3 = map (matches ex3) ["abracadabra","abracadaaaaabra"]
> ex4 = leibniz 'r' ex3
> test4 = map (matches ex4) ["abracadabra","abacadabra","abracaaaaadaba","abacadaba"]

fox_l a makes regular expressions accept any prefix ending just before an a. Combining fox_l and fox_r allows you to search for all elements of a language just between a chosen pair of symbols.

In ex5 we build a regexp for certain types of statement in a simple language and then pull out from it the part corresponding to substrings within parentheses. The undifferentiated part of ex5 would normally be written as "x = ([0-9]+[0-9]); *".

> ex5 = simplify $ fox_r '(' $ fox_l ')' $
>   let digits = Star (Sum (map Symbol "0123456789")) in
>        re "x = (" <*> digits <*> re "+" <*> digits <*> re ");" <+> Star (re " ")
> test5 = map (matches ex5) ["1+2","7-4"]

Open Question: There are more derivatives that can be constructed from the parts I've given. What do they do? What other multiplicative homomorphisms can you define? What derivatives do they give? (Hint: consider the property "R matches only a finite language").

Oh...I almost forgot the whole point of this. It's now possible to see that McBride's l-operator is a version of the Brzozowski derivative. In fact, the theorems he proves about l are pretty similar to the theorems that get proved about the Brzozowski derivative such as Theorem 4.4 here.

---

Some helper functions:

splits x returns all the ways of splitting x into pairs u and v such that u ++ v == x and v is non-empty.

> splits [] = []
> splits (a:as) = ([],a:as) : map f (splits as) where
>   f (x,y) = (a:x,y)

Some (incomplete) simplification rules for regular expressions

> simplify a = let b = simplify' a in if a==b then b else simplify b where
>   simplify' (Prod [a]) = a
>   simplify' (Prod a) | zero `elem` a = zero
>   simplify' (Prod (a:b)) = case a of
>       Prod x -> Prod $ x ++ map simplify b
>       a -> Prod $ filter (/= one) $ map simplify (a:b)
>
>   simplify' (Sum [a]) = a
>   simplify' (Sum (a:b)) = case a of
>       Sum x -> Sum $ x ++ map simplify b
>       a -> Sum $ nub $ sort $ filter (/= zero) $ map simplify (a:b)
>
>   simplify' (Star a) = Star (simplify a)
>   simplify' a = a

Fox's Ubiquitous Free Derivative

A post in three acts.

1. Regular Expressions

Every regular expression defines a language - the set of strings exactly matching the expression. Given a language L made from the symbols in S, define ∂_sL, where s is a symbol in S, to be the language made by taking all of the prefixes in L that end just before an s, deleting the final s. For example, if the language were L={MADAM} then ∂_ML={∅,MADA}, ∂_AL={M,MAD} and ∂_DL={MA}. We can write the regular expressions for these languages as 1+MADA, M+MAD and MA. So in this case, applying ∂_s trivially gives another regular language. So here's a problem: is the set of prefixes of a regular language L, ending at s, also regular?

Let's use ∂_s to define an operation on regular expressions that computes the regular expression for the prefixes.

Consider a regular expression of the form ab, where a and b are regular expresions. A prefix for the language ending at symbol s either ands in the 'a' part, corresponding to something matching ∂_sa, or it ends in the b part in which case it matches a∂_sb.

But so far I've just mentioned finite languages. What about the * operator. Well a prefix in a^* must consist of a sequence of strings matching a but end as a prefix of a. Ie. we need ∂_s(a^*)=a*∂_sa.

So it's not hard to see (I hope) that

∂_ss = 1
∂_ss' = 0 for symbol s'≠s
∂_s(a+b) = ∂_sa+∂_sb
∂_s(ab) = ∂_sa+a∂_sb
∂_s(a^*)=a*∂_sa

works nicely and proves that the answer to the above question is yes.

2. Dissecting types

Rather than regurgitate someone else's words, take a peek at Dissecting Data Structures. Define the operator, acting on types, ∂, by ∂ F X = Δ F X 1, where Δ is my poor man's version of del with a vertical line through it. From the definition of Δ you can see it satisfies:

∂(K A) X = 0
∂Id X = 1
∂(F+G) X=∂F X+∂G X
&part(F×G) X=∂F X×G 1+F X×∂G X

If you read McBride's paper, you'll see that the ∂ operator slices a type in half, throwing away the elements in the right half, keeping just the 'shape' of the right half.

3. Fox's Free Derivatives

Check out the definition of the Fox derivative in section 4 of Free Differential Calculus. The Fox derivative is used to compute the Alexander polynomial of a knot.

But I hope you can see that we have yet another definition that looks very similar. The common theme is differentiation with a twist. When you differentiate a product fg, one term is the usual one, fD(g), but the other term isn't D(f)g but the g is replaced by a drastically simplified version (either 1, or g(1)). Both versions 1 and 2 have a notion of chopping something in half at some point and keeping the left half. (Actually, you could tweak the definition in part 2 so that g is replaced by 1 too, then you'd no longer preserve the shape of the right half of the type and you'd really have a slicing operator, and the analogies would become clearer.)

I just need to relate the * operator in part 1 to the derivative of the inverse in the free derivative. Let me do this in a hand-waving manner: Conceptually, a^*=1+a+a²+… which might be considered the same as 1/(1-a). ∂(1/(1-a)) = (using 4(ii)) 1/(1-a)∂a = a^*∂a.

So what does it all mean? I've no idea. I'm just pointing out that the Fox derivative, or something like it, has what I think are hitherto unnoticed applications to computer science. Anyway, some of this will become clearer when I post some Haskell code implementing no less than five different types of derivatives of regular expressions (some time in the next week when I have a spare moment), and that should partially answer my question that I posed on λtU.

On further reflection, I can say something about what this means. The way to view derivatives is as an operation that chops things in 'half', destroying something in the middle. The Leibniz rules corresponds to the notion that when you chop a product you must either chop in the left side of the product or the right side. In the usual Leibniz rule you reassemble the pieces on either side - leaving a hole. With the Fox derivative you throw away what's on the right. With the ∂ operator in part 2 you keep only the shape on the right. And McBride's Δ operator gives you a choice of what do do with the left and right hand sides.

Oriented Fish and Hairy Balls

I have to admit, I have a tendency to write about abstract nonsense here. But I do so for entirely practical reasons - I don't want to run afoul of my employer by mentioning anything that could be construed as a trade secret. But today I thought I'd mention a small application of topology to computer graphics. The reason I can safely mention it is that I won't be revealing a technique or algorithm. Instead, I want to point out the opposite: a proof that a certain algorithm cannot possibly exist.

From time to time, anyone who works in computer graphics is likely to want to implement a function with this specification:

Implement a function f:R³→R³ such that for all non-zero x, f(x) is also non-zero and is perpendicular to x.

This kind of function is ubiquitous in graphics and most math libraries geared towards graphics that I have seen have some kind of implementation of this. A typical application is as follows: suppose you're animating shoals of fish. Using a variety of methods its easy to generate nice paths that guide the motion of thousands of fish. For example you could assume that the motion of the fish is approximately described by a divergence-free vector field (so that the fish don't 'pile up' anywhere) and then integrate these fields to give streamlines that the fish can follow. Unfortunately there is a slight catch with this. If you have some 3D geometry representing the fish's shape, knowing the position of the fish isn't enough to specify how it looks. You also need to give its orientation. An obvious choice is to orient the fish so that the line from tail to nose points along the flow. But you still need to be able to describe its orientation around this axis. To do this, we need to complete this orientation vector along the fish to an orthogonal basis ie. find a set of mutually orthogonal unit vectors, e₁, e₂, and e₃ such that e₁ points from tail to nose. We can easily define e₃ to be e₁×e₂, so if we had the function f above we could write e₂=f(e₁).

The function f appears in a multitude of other contexts too. It's useful when orienting 'virtual' cameras. For example the user might want to specify that the camera points at a target and have the graphics system select an orientation with this property.

However, anyone who's met these problems immediately hits an issue. Try writing code to do this for yourself if you haven't already. Chances are you'll end up with a conditional branch in your code, and when you flip from one branch to another you'll see it clearly causing ugly flips in the animation. So the above spec needs to be revised to:

Implement a continuous function f:R³→R³ such that for all non-zero x, f(x) is also non-zero and is perpendicular to x.

And here's the rub: the Hairy Ball Theorem says there is no such function. So if you're having trouble writing one, it's not because you're not being ingenious enough. It just can't be done! (Though there are lots of other solutions to the animation problems I mention above.) See, topology is useful after all, it can save you wasting your time on impossible tasks.

(Funnily enough though, the best way to see that this is true is using a suitable functor...)

Actually, I think anyone who works in graphics ought to read up on the Hairy Ball Theorem, the Ham Sandwich Theorem and the Borsuk-Ulam theorem as well as learning about the double cover of SO(3) (which I'll write about here in the near future).

That was a public service announcement.

Tying Knots Generically

A comment by a friend and a puzzle on Lambda the Ultimate made me realise that the function loeb that I defined a few weeks back does have a nice application.

> import Maybe
> import Data.Map

> loeb :: Functor a => a (a x -> x) -> a x
> loeb x = fmap (\a -> a (loeb x)) x

The problem it helps to solve is this: building circular datastructures in Haskell, otherwise known as Tying the Knot. loeb solves this problem in a very generic way that imposes the minimum of restrictions on the user.

Remember how in my previous article I showed how you could view loeb as a kind of spreadsheet evaluator. You have a container of functions and that is converted into another container by applying these functions to the final container. This is a circular definition, we're applying the functions to the final result to get the final result. If we think of the final container as a collection of answers to a bunch of problems, loeb allows each element to be constructed using references to any of the other answers in the container, even if they haven't actually been found yet. This is exactly what we need to make circular references.

Firstly, let's look at doing something similar to the Tying the Knot web page to solve Derek's puzzle. Here potential links are made by using a dictionary structure to allow names of objects to stand for direct references to objects. So define an operator to perform pointwise concatenation of string-valued functions:

> f +++ g = \x -> f x ++ g x

Now define lookup into the dictionary, substituting a default value if needed:

And now we can translate some example equations into Haskell by hand. I'm considering define a x *b and define b y *a. Note how it's trivial to write a parser to make these objects because there's nothing circular going on:

> equations = fromList
>     [
>         ("a",const "x"+++lookup "b"),
>         ("b",const "y"+++lookup "a")
>     ] where lookup = Data.Map.findWithDefault ""

We now just apply loeb to get a dictionary containing the actual values:

> values = loeb equations
> a = fromJust $ Data.Map.lookup "a" values

We can use a similar method to make a doubly linked list. Here is "A", "B" and "C" in a list with each letter pointing to its predecessor and successor if it exists:

> data DList a = Nil | Node a (DList a) (DList a) deriving Show
> 
> val (Node a _ _) = a
> rt (Node _ _ r) = r
> lft (Node _ l _) = l
> 
> nodes = fromList
>    [
>        ("a", \d -> Node "A" Nil (lookup "b" d)),
>        ("b", \d -> Node "B" (lookup "a" d) (lookup "c" d)),
>        ("c", \d -> Node "C" (lookup "b" d) Nil)
>    ] where lookup = Data.Map.findWithDefault Nil
> 
> abc = loeb nodes
> b = fromJust $ Data.Map.lookup "b" abc

b is now a doubly linked list and we can use lft and rt to walk up and down it.

And here's a nice advantage of using loeb. You don't have to restrict yourself to storing information about links in a dictionary. Any functorial container will do. For example here's a doubly linked circular list of the integers from 0 to 99 made without using any names:

> zlist = [\z -> Node i (z!!((i-1) `mod` 100)) (z!!((i+1) `mod` 100)) | i <- [0..99]]
> z = loeb zlist

Note how working your way up and down zlist is a slow process because you need to use the slow (!!) to access elements. But stepping up and down any element of z takes time O(1). Potentially very useful!

But there is a slight problem with loeb. I'd like to be able to handle errors better and the obvious way to do that is by allowing the use of monads to handle exceptions, for example if a name lookup in a dictionary fails. I attempted to do this as follows:

> class Monad m => MFunctor m f where
>     fmapM :: (a -> m b) -> f a -> m (f b)

> instance Monad m => MFunctor m [] where
>     fmapM = mapM

> loebM :: (Monad m,MFunctor m a) => a (a x -> m x) -> m (a x)
> loebM x = r where
>     r = do
>         lx <- r
>         fmapM (\a -> a lx) x

> u = loebM [const (Just 1),const (Just 2)]

But attempting to evaluate u overfows the stack. It also failed when I tried modifying it to use MonadFix. So here's a puzzle. Can you write a monadic version of loeb? It seems like something that would be very useful for writing interpreters for programming languages that allow mutually recursive definitions.

So my original intuition a few weeks back was right. Taking an interesting axiom of modal logic, interpreting as a type, and them finding an implementation, did lead to something useful.

---

By the way, while I was thinking about loeb I was led to consider what I call tautological containers. These are containers where each element tells you how to extricate it from the container. For example, these can be thought of as tautological: (fst,snd), [head,head.tail], [(!!n) | n <- [0..]] and [("a",fromJust . lookup "a"),("b",fromJust . lookup "b")]. Can you see any use for these objects? Are they universal for some category diagram? I have partial answers to these questions but maybe you can make more sense of them than me.

---

Update: Ocaml version by Matias Giovannini.

Evaluating cellular automata is comonadic

Paul Potts's post inspired me to say something about cellular automata too.

So here's the deal: whenever you see large datastructures pieced together from lots of small but similar computations there's a good chance that we're dealing with a comonad. In cellular automata we compute the value of each cell in the next generation by performing a local computation based on the neighbourhood of that cell. So cellular automata look like they might form a good candidate for comonadic evaluation.

I want to work on 'universes' that extend to infinity in both directions. And I want this universe to be constructed lazily on demand. One way of doing that is to represent a 'universe' as a centre point, a list of all elements to the left of that centre point and a list of all elements to the right. Here's a suitable
type:

> data U x = U [x] x [x]

For example U [-1,-2..] 0 [1,2..] can be thought of as representing all of the integers in sequence.

But this actually contains slightly more information than a list that extends to infinity both ways. The centre point forms a kind of focus of attention. We could shift that focus of attention left or right. For example consider

U [-2,-3..] (-1) [0,1..]

This represents the same sequence of integers but the focus has been shifted left. So think of the type U x as being a doubly infinite sequence with a cursor. (In fact, this makes it a kind of zipper.)

We can formalise the notion of shifting left and right as follows:

> right (U a b (c:cs)) = U (b:a) c cs
> left  (U (a:as) b c) = U as a (b:c)

An object of type U is semantically like a C pointer into a const block of memory. You can increment it, decrement it and dereference it using the function I'll call coreturn below.

As U is a kind of list structure, it needs a map. In fact, we can define fmap for it:

> instance Functor U where
>    fmap f (U a b c) = U (map f a) (f b) (map f c)

Now the fun starts. First I'll bemoan the fact that Comonads aren't in the standard Haskell libraries (at least I don't think they are). So I have to define them myself:

> class Functor w => Comonad w where
>    (=>>)    :: w a -> (w a -> b) -> w b
>    coreturn :: w a -> a
>    cojoin     :: w a -> w (w a)
>    x =>> f = fmap f (cojoin x)

cojoin is the dual to the usual join function. I've chosen to do things the category theoretical way and define =>> in terms of cojoin.

And here's why U forms a Comonad:

> instance Comonad U where
>    cojoin a = U (tail $ iterate left a) a (tail $ iterate right a)
>    coreturn (U _ b _) = b

Look closely at cojoin. It turns a into a 'universe' of 'universes' where each element is a copy of a shifted left or right a number of times. This is where all the work is happening. The reason we want to do this is as follows: we want to write rules that work on the local neighbourhoods of our universe. We can think of a universe with the cursor pointing at a particular element as being an element with a neighbourhood on each side. For example, we can write a cellular automaton rule like this:

> rule (U (a:_) b (c:_)) = not (a && b && not c || (a==b))

In order to apply this everywhere in the universe we need to apply the rule to each possible shift of the universe. And that's what cojoin does, it constructs a universe of all possible shifts of a. Compare with what I said here. So believe it or not, we've already written the code to evaluate cellular automata. u =>> rule applies the rule to u. The rest is just boring IO:

> shift i u = (iterate (if i<0 then left else right) u) !! abs i
>
> toList i j u = take (j-i) $ half $ shift i u where
>    half (U _ b c) = [b] ++ c
>
> test = let u = U (repeat False) True (repeat False)
>       in putStr $
>          unlines $
>          take 20 $
>          map (map (\x -> if x then '#' else ' ') . toList (-20) 20) $
>          iterate (=>> rule) u

Lazy infinite structures, comonads, zippers. I think I'm just beginning to get the hang of this functional programming lark! Over Xmas I might feel ready to try writing a piece of code longer than a dozen or so lines.

Anyway, I must end with a credit. I probably wouldn't have come up with this if I hadn't read this paper by Uustalu and Vene.

Why isn't ListT [] a monad?

Well I'm back from my vacation. But this isn't my personal blog so I think it's time to dive right in. Anything to take my mind off these mosquito bites...

So consider the free semiring R generated by some set S. In other words, S consists of finite sums and products of 0, 1 and elements of S and where the only simplification rules allowed are

a+0=a, a+b=b+a, a+(b+c)=(a+b)+c

a1=1, a(bc) = (ab)c

a(b+c) = ab+ac and (a+b)c = ac+bc.

For example, consider a term like ab+c(ef+gh). We can use distributivity to multiply this out to ab+cef+cgh. There's not much else we can do to simplify this. Now notice that if you multiply out all of the brackets wherever you can, ie. use distributivity until you can no longer, you end up with an expression that is a sum of terms and each term is a product of generators from S. (The sum of terms may be empty, in which case it's zero, and each of the terms may be a product of zero generators, in which case it is 1). This allows us to write elements of R in a canonical form: as a set of terms where each term is an ordered list of generators. For example, we could write ab+cef+cgh as {[a,b],[c,e,f],[c,g,h]}. Notice how the commutativity of addition is represented by the fact that we're using a set of terms, but each term is a list of generators because we're making no assumption about the commutativity of multiplication.

In Haskell it's more convenient to work with lists. Se we'll represent our running example as [['a','b'],['c','e','f'],['c','g','h']]. So if S is the Haskell type corresponding to our set of generators, then [[S]] can be thought of as the free semiring generated by elements of S, with the proviso that we consider two elements equal if one can be reordered to the other. Note that []=0 and [[]] = 1.

Now suppose that S is itself a free semiring generated by T. Then R = [[ [[T]] ]], modulo ordering. If you think about it, there's a nice way to use the algebraic structure to 'flatten' an element of [[ [[T]] ]] down to an element of [[T]]. R is freely generated by elements of S, in other words it consists of sums of products of elements of S. But the elements of S are themselves sums and products. So we have sums of products of sums of products. I emphasised two words in that sentence. This is because R contains subparts that are products of sums. If we multiply these out in the usual way, we get sums of sums of products of products. Group the sums and products together, and we get back to sums of products. Here's an example: any element of T trivially gives an element of S. So if a,b,c,e,f,g and h are all in T, then a,b,c and ef+gh are all in S and hence are generators of R, so ab+c(ef+gh) is in R. Multiply out and we obviously have the element ab+cef+cgh of S. It's not hard to see how this generalises to map any element of S back down to T. So we have maps

T -> [[T]] (by trivial inclusion of generators)

[[ [[T]] ]] -> [[T]] (by above argument)

Look to you like a monad? Of course it does! :-) But what monad is it?

First, let's implement the map [[ [[T]] ]] -> [[T]]. After tinkering I came up with

> import Control.Monad.List

> flatten x = concatMap (map concat . sequence) x

We can test it out

> x = [[ [['a']],[['b']] ], [ [['c']],[['e','f'],['g','h']] ]]
> test1 = flatten x

It multiplies out exactly the way we want.

Now compare with running this:

> x' = ListT [[ ListT [['a']],ListT [['b']] ],
>           [ ListT [['c']],ListT [['e','f'],['g','h']] ]]
> test2 = runListT $ join x'

In other words, apart from the ListT fluff, join for ListT [] is flatten. So if we consider lists as a mechanism to represent sets, ListT [] is revealed as the monad of semirings. I'm not sure, but I think that historically this is where monads originally came from. Certainly there are many papers on the relationship between monads and algebraic structures like semirings.

And now I can answer my original question. In a semiring, addition is commutative, so the order of terms doesn't matter. But in ListT [], we're using lists, and order does matter in a list. So if we do take order into account, then really ListT [] is the monad of semirings where both addition and multiplication is non-commutative. And here's the problem: in general, there is no such thing as a freely generated semiring with non-commutative addition.

Here's why: consider the expression ((a+b)+(a+b))*((a+b)+(a+b)). Multiply out the inner parentheses first and we get

(a+b+a+b)*(a+b+a+b)

= a*a+a*b+a*a+a*b+…

Now multiply out the outer parentheses first and we get

(a+b)*(a+b)+(a+b)*(a+b)+(a+b)*(a+b)+(a+b)*(a+b)

= a*a+a*b+b*a+b*b+…

The terms are coming out in a different order. Essentially distributivity doesn't work in a semiring with non-commutative addition.

This translates directly into failure of one of the monad laws. First write our expression as an object of type ListT []:

u = a+b

> u = ListT [["a"],["b"]]

v = (a+b)+(a+b)

> v = ListT [[u],[u]]

w = ((a+b)+(a+b))*((a+b)*(a+b))

> w = ListT [[v,v]]

join multiplies out parentheses. So working from the outer parentheses inwards we can use:

> expanded1 = join $ join w
> go1 = runListT expanded1

Working outwards from the inner parentheses we get:

> expanded2 = join $ fmap join w
> go2 = runListT expanded2

(Note how I use fmap to 'duck down' through the outer layer of parentheses so as to multiply out each of the subexpressions first.)

Evaluate go1 and go2 and you'll see that they corresponds to the two ways of multiplying out that I gave above. And more importantly, the values of expanded1 and expanded2 aren't equal, meaning that join . join = join . fmap join isn't satisfied. You may recognise this: it's one of the monad laws. (At least it's one of the monad laws written the way category theorists write them.) So we ain't got no monad.

---

I think this is now a fairly complete analysis of what ListT [] is all about. So one obvious remaining question is: where do games come into all this? The answer is that games form a semiring in a way that I haven't seen documented anywhere (though is surely common knowledge). I'd explain but I've run out of time...

Note that the reason ListT [] isn't a monad is that [] isn't commutative, in some sense. This has been observed many times in the past. Two papers mentioning this are this and this.

I actually figured out all of this stuff before this. I realised that the trees I was scribbling on the backs of my envelopes to represent elements of semirings could actually be generalised to just about any kind of tree, so I wrote about the general case first.

I prefer my example of ListT [] failing to be a monad to the examples given here. The latter make use of the IO monad so they aren't quite as 'pure'.