Geometric Algebra for Free!
I'm a fan of Haskell. But I'm not a zealot. I'll be the first to tell you that some of the stuff down below works better in C++. But C++ has syntax uglier than a mandrill's rear end so I'll sacrifice functionality for readability. This is literate Haskell by the way. Copy and paste this whole entry into a .lhs file and run it with hugs -98 or ghci -fglasgow-exts.
The first bit of functionality I'll sacrfice is some of the genericity of my code so let's start with fixing a base field. Everything I say below will be over the reals:
Anyway, today's theme is about how writing your code generically allows you to combine classes together in unanticipated ways to give you useful stuff - more or less for free.
If you've done some mathematical coding you probably have a complex number class lying around somewhere. If it's written in a modern language it's probably generic allowing you to choose your base structure. For example in C++ you might be able to use complex<float> or complex<int>. If you don't have such a class you can have a Haskell one for free:
The usual complex numbers are given by Complex R. But here's a question: does Complex (Complex R) mean anything. Complex a requires a to be a Num, but Complex R is a Num. So Complex (Complex R) compiles fine. But does it mean anything mathematically? It does, it corresponds to C⊗C. (Remember, I'm working over the reals, so by ⊗ I mean ⊗R.) If you're shaky on tensor products then I can give a one line summary for the case of working over the reals: A⊗B means opening up the object of type A and replacing every occurence of a real number in it with a copy of a B. For example, a complex number can be thought of as an object that contains two real numbers. So C⊗C is an object that contains two complex numbers and which multiplies in a certain way. So the type constructor Complex actually corresponds quite well with the functor _⊗C. It's even better, if ⊗ is thought of as the tensor product in the category of algebras, not just vector spaces, then the Complex type constructor computes a tensor product of algebras. This is kinda obvious but I've never seen anyone say this explicitly. (Clearly the same is true of a C++ template like complex<>.)
So let's look at some examples. Maybe you also have a matrix class lying around. If you don't, here's a 2 by 2 matrix class for you:
We can start putting our words into practice. Use the notation A(n) to mean the n by n matrices with entries in A. Then A(n)=R(n)⊗A. So the type constructor Matrix is the functor _(2). Another theorem is that R(m)⊗R(n)=R(mn). This means that Matrix (Matrix R) gives 4 by 4 matrices. Neat eh? We're already getting some classes for free. But we can get bigger stuff for free.
I'm also a fan of geometric algebra. But I'm not a zealot. Geometric algebra can do a lot of neat stuff but I'm not one of those people campaigning to get everyone to throw out their old vectors and tensors to replace them with geometric algebra.
Geometric algebra is all about Clifford algebras. The Clifford algebra, Cliff(n), is the real numbers extended by throwing in n elements ei, for i from 1 to n, and all of the things you can get my multiplying and adding these numbers. In addition we impose the rules that ei²=-1 and eiej=-ejei. You already know some Clifford algebras. Cliff(0) is just r. Cliff(1) is C because i serves perfectly well as e1. Cliff(2) is probably familiar to you as well, it's just the quaternions. Ah...but quaternions have three values whose square is -1, so surely it's Cliff(3). Well, no actually. e1=i, e2=j and k is reallly just ij, ie. e1e2, and so isn't something new. In general, Cliff(n) is 2n dimensional with a basis consisting of all of the possible products of the ei. Anyway, if you don't already have a Quaternion class, here's one for you:
There are also another bunch of Clifford algebras called Cliff'(n). Cliff'(n), is the real numbers extended by throwing in n elements ei, for i from 1 to n. With the properties ei²=1 and eiej=-ejei. Note the 1 instead of -1. Cliff'(1) is known as the split complex numbers. You might not have an implementation of these, so here you go, gratis:
You already have an implementation of Cliff'(2). Believe it or not, it's isomorphic to the 2 by 2 matrices. e1=[[-1,0],[0,1]] and e2=[[0,1],[1,0]], using the obvious Haskell lists to represent matrices. (Try squaring those matrices and seeing what you get.)
That's all well and good, but what about Cliff(n) for larger n. The obvious thing is to represent the elements as arrays of length 2n. But to multiply these things we'd have to loop over each element of the first multiplicand and then over each element of the second meaning a loop that takes time 22n. But who wants to write code when you can get stuff for free? Check outh this paper on the classification of Clifford algebras. On page 5 there's an important theorem
Cliff(n+2)=Cliff(2)⊗Cliff'(n)
Cliff'(n+2)=Cliff'(2)⊗Cliff(n)
Now we know how to implement tensor products this becomes trivial. In Haskell notation,
and so on.
That's it. We've built Cliff(n) and Cliff'(n) out of a bunch of simple classes that you might already have had lying around. Geometric algebra for free! It almost seems sad to see a beautiful classification theorem demystified in such a simple way. Still, you don't see pairs of mutually recursive towers of types like this every day.
But actually this is all a con. It's no good having something isomorphic to what we want without the actual isomorphism. For example where is e3 in Cliff4', say? So we now have to do a tiny bit of real work. Start with a type class:
i is the function that embeds our base field, the reals. And e i a is intended to be aei.
And now we simply assert:
Now we can define i and e for all of the algebras. For the details look at the proof of Theorem 3.4. They translate almost directly as:
And now we really do have usable Clifford algebras. Here are some more defintions:
(In C++ you can define 'metafunctions' that map integers to classes. Doing something similar in Haskell is pretty ugly.)
And I've learnt my lesson from my monad tutorial. Here's a little test to make sure I'm actually telling the truth:
If fulltest evaluates to True then we know that for Cliff9 and Cliff9' we have created Haskell objects that satisfy the defining properties for ei.
The first bit of functionality I'll sacrfice is some of the genericity of my code so let's start with fixing a base field. Everything I say below will be over the reals:
> type R = Float
Anyway, today's theme is about how writing your code generically allows you to combine classes together in unanticipated ways to give you useful stuff - more or less for free.
If you've done some mathematical coding you probably have a complex number class lying around somewhere. If it's written in a modern language it's probably generic allowing you to choose your base structure. For example in C++ you might be able to use complex<float> or complex<int>. If you don't have such a class you can have a Haskell one for free:
> data Complex a = C a a deriving (Eq,Show)
> instance Num a => Num (Complex a) where
> fromInteger a = C (fromInteger a) 0
> (C a b)+(C a' b') = C (a+a') (b+b')
> (C a b)-(C a' b') = C (a-a') (b-b')
> (C a b)*(C a' b') = C (a*a'-b*b') (a*b'+b*a')
The usual complex numbers are given by Complex R. But here's a question: does Complex (Complex R) mean anything. Complex a requires a to be a Num, but Complex R is a Num. So Complex (Complex R) compiles fine. But does it mean anything mathematically? It does, it corresponds to C⊗C. (Remember, I'm working over the reals, so by ⊗ I mean ⊗R.) If you're shaky on tensor products then I can give a one line summary for the case of working over the reals: A⊗B means opening up the object of type A and replacing every occurence of a real number in it with a copy of a B. For example, a complex number can be thought of as an object that contains two real numbers. So C⊗C is an object that contains two complex numbers and which multiplies in a certain way. So the type constructor Complex actually corresponds quite well with the functor _⊗C. It's even better, if ⊗ is thought of as the tensor product in the category of algebras, not just vector spaces, then the Complex type constructor computes a tensor product of algebras. This is kinda obvious but I've never seen anyone say this explicitly. (Clearly the same is true of a C++ template like complex<>.)
So let's look at some examples. Maybe you also have a matrix class lying around. If you don't, here's a 2 by 2 matrix class for you:
> data Matrix a = M a a a a deriving (Eq,Show)
> instance Num a => Num (Matrix a) where
> fromInteger a = M (fromInteger a) 0 0 (fromInteger a)
> (M a b c d)+(M a' b' c' d') = M (a+a') (b+b')
> (c+c') (d+d')
> (M a b c d)-(M a' b' c' d') = M (a-a') (b-b')
> (c-c') (d-d')
> (M a b c d)*(M a' b' c' d') = M (a*a'+b*c') (a*b'+b*d')
> (c*a'+d*c') (c*b'+d*d')
We can start putting our words into practice. Use the notation A(n) to mean the n by n matrices with entries in A. Then A(n)=R(n)⊗A. So the type constructor Matrix is the functor _(2). Another theorem is that R(m)⊗R(n)=R(mn). This means that Matrix (Matrix R) gives 4 by 4 matrices. Neat eh? We're already getting some classes for free. But we can get bigger stuff for free.
I'm also a fan of geometric algebra. But I'm not a zealot. Geometric algebra can do a lot of neat stuff but I'm not one of those people campaigning to get everyone to throw out their old vectors and tensors to replace them with geometric algebra.
Geometric algebra is all about Clifford algebras. The Clifford algebra, Cliff(n), is the real numbers extended by throwing in n elements ei, for i from 1 to n, and all of the things you can get my multiplying and adding these numbers. In addition we impose the rules that ei²=-1 and eiej=-ejei. You already know some Clifford algebras. Cliff(0) is just r. Cliff(1) is C because i serves perfectly well as e1. Cliff(2) is probably familiar to you as well, it's just the quaternions. Ah...but quaternions have three values whose square is -1, so surely it's Cliff(3). Well, no actually. e1=i, e2=j and k is reallly just ij, ie. e1e2, and so isn't something new. In general, Cliff(n) is 2n dimensional with a basis consisting of all of the possible products of the ei. Anyway, if you don't already have a Quaternion class, here's one for you:
> data Quaternion a = Q a a a a deriving (Eq,Show)
> instance Num a => Num (Quaternion a) where
> fromInteger a = Q (fromInteger a) 0 0 0
> (Q a b c d)+(Q a' b' c' d') = Q (a+a') (b+b')
> (c+c') (d+d')
> (Q a b c d)-(Q a' b' c' d') = Q (a-a') (b-b')
> (c-c') (d-d')
> (Q a b c d)*(Q a' b' c' d') = Q (a*a'-b*b'-c*c'-d*d')
> (a*b'+b*a'+c*d'-d*c')
> (a*c'+c*a'+d*b'-b*d')
> (a*d'+d*a'+b*c'-c*b')
There are also another bunch of Clifford algebras called Cliff'(n). Cliff'(n), is the real numbers extended by throwing in n elements ei, for i from 1 to n. With the properties ei²=1 and eiej=-ejei. Note the 1 instead of -1. Cliff'(1) is known as the split complex numbers. You might not have an implementation of these, so here you go, gratis:
> data Split a = S a a deriving (Eq,Show)
> instance Num a => Num (Split a) where
> fromInteger a = S (fromInteger a) (fromInteger a)
> (S a b)+(S a' b') = S (a+a') (b+b')
> (S a b)-(S a' b') = S (a-a') (b-b')
> (S a b)*(S a' b') = S (a*a') (b*b')
You already have an implementation of Cliff'(2). Believe it or not, it's isomorphic to the 2 by 2 matrices. e1=[[-1,0],[0,1]] and e2=[[0,1],[1,0]], using the obvious Haskell lists to represent matrices. (Try squaring those matrices and seeing what you get.)
That's all well and good, but what about Cliff(n) for larger n. The obvious thing is to represent the elements as arrays of length 2n. But to multiply these things we'd have to loop over each element of the first multiplicand and then over each element of the second meaning a loop that takes time 22n. But who wants to write code when you can get stuff for free? Check outh this paper on the classification of Clifford algebras. On page 5 there's an important theorem
Now we know how to implement tensor products this becomes trivial. In Haskell notation,
> type Cliff1 = Complex R
> type Cliff1' = Split R
> type Cliff2 = Quaternion R
> type Cliff2' = Matrix R
> type Cliff3 = Quaternion Cliff1'
> type Cliff3' = Matrix Cliff1
> type Cliff4 = Quaternion Cliff2'
> type Cliff4' = Matrix Cliff2
and so on.
That's it. We've built Cliff(n) and Cliff'(n) out of a bunch of simple classes that you might already have had lying around. Geometric algebra for free! It almost seems sad to see a beautiful classification theorem demystified in such a simple way. Still, you don't see pairs of mutually recursive towers of types like this every day.
But actually this is all a con. It's no good having something isomorphic to what we want without the actual isomorphism. For example where is e3 in Cliff4', say? So we now have to do a tiny bit of real work. Start with a type class:
> class Clifford b where
> i :: R -> b
> e :: Int -> R -> b
i is the function that embeds our base field, the reals. And e i a is intended to be aei.
And now we simply assert:
> instance Clifford R where
> i a = a
> e _ _ = error ""
> instance Clifford (Complex R) where
> i a = C a 0
> e 1 a = C 0 a
> instance Clifford (Split R) where
> i a = S a a
> e 1 a = S a (-a)
Now we can define i and e for all of the algebras. For the details look at the proof of Theorem 3.4. They translate almost directly as:
> instance (Num b,Clifford b) => Clifford (Quaternion b) where
> i a = Q (i a) 0 0 0
> e 1 a = Q 0 (i a) 0 0
> e 2 a = Q 0 0 (i a) 0
> e n a = Q 0 0 0 (e (n-2) a)
> instance (Num b,Clifford b) => Clifford (Matrix b) where
> i a = let b = i a in M b 0 0 b
> e 1 a = let b = i a in M (-b) 0 0 b
> e 2 a = let b = i a in M 0 b b 0
> e n a = let b = e (n-2) a in M 0 b (-b) 0
And now we really do have usable Clifford algebras. Here are some more defintions:
> type Cliff5 = Quaternion Cliff3'
> type Cliff5' = Matrix Cliff3
> type Cliff6 = Quaternion Cliff4'
> type Cliff6' = Matrix Cliff4
> type Cliff7 = Quaternion Cliff5'
> type Cliff7' = Matrix Cliff5
> type Cliff8 = Quaternion Cliff6'
> type Cliff8' = Matrix Cliff6
> type Cliff9 = Quaternion Cliff7'
> type Cliff9' = Matrix Cliff7
(In C++ you can define 'metafunctions' that map integers to classes. Doing something similar in Haskell is pretty ugly.)
And I've learnt my lesson from my monad tutorial. Here's a little test to make sure I'm actually telling the truth:
> t i = e i 1 :: Cliff9
> t' i = e i 1 :: Cliff9'
>
> test a t n = map (\(i,j) -> (t i)*(t j)+(t j)*(t i)==0) [(i,j) | i <- [1..n], j <- [1..(i-1)]] ++
> map (\i -> (t i)*(t i)==a) [1..n]
> fulltest = and $ (test (-1) t 9) ++ (test 1 t' 9)
If fulltest evaluates to True then we know that for Cliff9 and Cliff9' we have created Haskell objects that satisfy the defining properties for ei.
One last thing, let's consider the complexity of this code. The naive algorithm multiplies two elements of Cliff(n) in time O(22n). It's an easy exercise to see that the multiplication time is O(27n/4). Not bad for something that was free. We can acualy do quite a bit better with some small changes - for example by using Strassen-like tricks for all but the Split multiplications.
But if you really want to do this stuff properly check out this paper which I haven't yet read properly.
I'd also like to add that Cliff(n+8)=R(16)⊗Cliff(n), so we have a kind of repeat with period 8. This is known as Bott periodicity and it lies at the heart of many of the most beautiful things in mathematics including some fascinating properties of the number 8.
Labels: haskell, mathematics
posted by sigfpe at Wednesday, August 30, 2006
24 Comments:
It would be great to have your blog on planet.haskell.org . You can email antti [dash] juhani [at] kaijanaho.info to confirm if you like, as he won't add you on without your permission.
Dude, you rock. Very interesting stuff.
When you say some things are better described in C++, I presume you mean the function that takes an n and generates the type Cliff(n).
But just to show that how it looks in Haskell here's some code. The big difference between the Haskell and C++ versions of doing this is that the Haskell version is actually type checked, whereas in C++ it would only be checked when instantiating the template. (Also the Haskell version allows the argument n to be a dynamic (runtime) argument rather than static.)
First some preliminaries no define type level natural numbers.
> data Zero
> data Succ a
> type N0 = Zero
> type N1 = Succ N0
> type N2 = Succ N1
> type N3 = Succ N2
> type N4 = Succ N3
> type N5 = Succ N4
> type N6 = Succ N5
> type N7 = Succ N6
> type N8 = Succ N7
> type N9 = Succ N8
> type Plus2 n = Succ (Succ n)
The MkCliff and MkCliff' classes are really functions from n, a type level nutural, to c, which is Cliff(n) resp. Cliff'(n).
> class MkCliff n c | n -> c
> class MkCliff' n c | n -> c
The base cases for the Clifford algebras.
> instance MkCliff N0 R
> instance MkCliff' N0 R
> instance MkCliff N1 (Complex R)
> instance MkCliff' N1 (Split R)
And finally the inductive cases.
> instance (MkCliff' n c) => MkCliff (Plus2 n) (Quaternion c)
> instance (MkCliff n c) => MkCliff' (Plus2 n) (Matrix c)
Define the i and e functions, but now taking an extra argument that define which algebra to embed into.
> mk_i :: (MkCliff n c, Clifford c) => n -> R -> c
> mk_i _ = i
> mk_e :: (MkCliff n c, Clifford c) => n -> Int -> R -> c
> mk_e _ = e
> mk'_i :: (MkCliff' n c, Clifford c) => n -> R -> c
> mk'_i _ = i
> mk'_e :: (MkCliff' n c, Clifford c) => n -> Int -> R -> c
> mk'_e _ = e
And then define t and t' again.
> n9 = undefined :: N9
> t i = mk_e n9 i 1
> t' i = mk'_e n9 i 1
augustss,
That's pretty nice. I particularly liked the use of 'undefined' that almost makes the type N9 feel like a first class object in Haskell. I'd not seen that trick before.
One thing I couldn't figure out was how to make all of the code generic with respect to the base field - so I locked it down with 'type R = Float'. I tried lots of messing about with multiparameter type classes and fundeps but a compile error would always crop up somewhere.
And I'm kinda surprised that someone actually did something with my code!
There's more than one way to skin a cat, but here's one way to make the base type a parameter.
I've also introduced a new method in Clifford that returns the grade so the test is more generic.
The use of I and Compose is a bit of an ugly trick since Haskell only allows data types in instance declarations. There must be a neater way to do this.
It's interesting how MkCliff only builds the data type structure of the algebra and then the Clifford class provides the computational contents.
class Clifford b r where
i :: r -> b r
e :: Int -> r -> b r
g :: b r -> Int -- the grade
newtype I a = I a
instance Clifford I a where
i a = I a
e _ _ = error ""
g _ = 0
instance (Num a) => Clifford Complex a where
i a = C a 0
e 1 a = C 0 a
g _ = 1
instance (Num a) => Clifford Split a where
i a = S a a
e 1 a = S a (-a)
g _ = 1
newtype Compose f g a = O (f (g a)) deriving (Eq, Show, Num)
instance (Num (b r), Clifford b r) => Clifford (Compose Quaternion b) r where
i a = O $ Q (i a) 0 0 0
e 1 a = O $ Q 0 (i a) 0 0
e 2 a = O $ Q 0 0 (i a) 0
e n a = O $ Q 0 0 0 (e (n-2) a)
g ~(O (Q a _ _ _)) = 2 + g a
instance (Num (b r), Clifford b r) => Clifford (Compose Matrix b) r where
i a = O $ let b = i a in M b 0 0 b
e 1 a = O $ let b = i a in M (-b) 0 0 b
e 2 a = O $ let b = i a in M 0 b b 0
e n a = O $ let b = e (n-2) a in M 0 b (-b) 0
g ~(O (M a _ _ _)) = 2 + g a
class MkCliff n (c :: * -> *) | n -> c
class MkCliff' n (c :: * -> *) | n -> c
instance MkCliff N0 I
instance MkCliff' N0 I
instance MkCliff N1 Complex
instance MkCliff' N1 Split
instance (MkCliff' n c) => MkCliff (Plus2 n) (Compose Quaternion c)
instance (MkCliff n c) => MkCliff' (Plus2 n) (Compose Matrix c)
mk_i :: (MkCliff n c, Clifford c a) => n -> a -> c a
mk_i _ = i
mk_e :: (MkCliff n c, Clifford c a) => n -> Int -> a -> c a
mk_e _ = e
mk'_i :: (MkCliff' n c, Clifford c a) => n -> a -> c a
mk'_i _ = i
mk'_e :: (MkCliff' n c, Clifford c a) => n -> Int -> a -> c a
mk'_e _ = e
n9 = undefined :: N9
t i = mk_e n9 i (1::Double)
t' i = mk'_e n9 i (1::Double)
test a t = map (\(i,j) -> (t i)*(t j)+(t j)*(t i)==0) [(i,j) | i <- [1..n], j <- [1..(i-1)]] ++
map (\i -> (t i)*(t i)==a) [1..n]
where n = g (t undefined)
fulltest = and $ (test (-1) t) ++ (test 1 t')
Hey! That was quick work!
Now you're using Haskell features I'm not familiar with, like the line "class MkCliff n (c :: * -> *) | n -> c". It's going to take me some time to digest this.
I think Clifford algebras might make a great case study for developing ideas about statically typed programming languages specifically for mathematics. Haskell is good, but you're having to use ingenious tricks to express things that mathematicians would like to be easy.
The (c :: * -> *) is just there to tell the type checker that the type variable c has kind *->*. The class declaration itself does not contain enough information to determine the kind otherwise, so the compiler would assume *.
The 'n -> c' part just says that n uniquely determines c.
Yes, many things that are simple in mathematics get more complicated in programming languages. The big reason being that mathematicians are smarter than computers. :)
A few more comments:
First,
Why is 'i' a separate method? Wouldn't 'e 0' instead of 'i' make sense?
Second,
The code you wrote can be parametrized on R without any tricks at all, see below. The only troublesome case is 'instance Clifford R'. When you make R a parameter this gets to be a very general instance declaration, and those are often "dangerous", by that I mean that they apply far to often. So I've introduced the I type that is just the identity.
Furthermore, the 'b -> r' part of the class declaration is only important if you want to avoid adding a lot of type signatures. It just says that the base type, r, is uniquely determined by the Clifford type itself (which is true). Without this vital clue the compiler sometimes cannot deduce types.
> type Cliff1 r = Complex r
> type Cliff1' r = Split r
> type Cliff2 r = Quaternion r
> type Cliff2' r = Matrix r
> type Cliff3 r = Quaternion (Cliff1' r)
> type Cliff3' r = Matrix (Cliff1 r)
> type Cliff4 r = Quaternion (Cliff2' r)
> type Cliff4' r = Matrix (Cliff2 r)
> class Clifford b r | b -> r where
> i :: r -> b
> e :: Int -> r -> b
> newtype I r = I r deriving (Eq, Ord, Show, Num)
> instance Clifford (I r) r where
> i a = I a
> e _ _ = error ""
> instance (Num r) => Clifford (Complex r) r where
> i a = C a 0
> e 1 a = C 0 a
> instance (Num r) => Clifford (Split r) r where
> i a = S a a
> e 1 a = S a (-a)
> instance (Num b,Clifford b r) => Clifford (Quaternion b) r where
> i a = Q (i a) 0 0 0
> e 1 a = Q 0 (i a) 0 0
> e 2 a = Q 0 0 (i a) 0
> e n a = Q 0 0 0 (e (n-2) a)
> instance (Num b,Clifford b r) => Clifford (Matrix b) r where
> i a = let b = i a in M b 0 0 b
> e 1 a = let b = i a in M (-b) 0 0 b
> e 2 a = let b = i a in M 0 b b 0
> e n a = let b = e (n-2) a in M 0 b (-b) 0
> type Cliff5 r = Quaternion (Cliff3' r)
> type Cliff5' r = Matrix (Cliff3 r)
> type Cliff6 r = Quaternion (Cliff4' r)
> type Cliff6' r = Matrix (Cliff4 r)
> type Cliff7 r = Quaternion (Cliff5' r)
> type Cliff7' r = Matrix (Cliff5 r)
> type Cliff8 r = Quaternion (Cliff6' r)
> type Cliff8' r = Matrix (Cliff6 r)
> type Cliff9 r = Quaternion (Cliff7' r)
> type Cliff9' r = Matrix (Cliff7 r)
> t i = e i 1 :: (Cliff9 Float)
> t' i = e i 1 :: (Cliff9' Float)
augustss,
First,
Why is 'i' a separate method? Wouldn't 'e 0' instead of 'i' make sense?
I definitely want to keep i and e separate. Conceptually they do different things: one is embedding the base field in the algebra and the other is embedding an n-dimensional vector space in the algebra. In a suitable programming language I'd replace e with something of type e :: Vector(n) -> Cliff(n). I definitely don't want to fold in the scalar part to get Vector(n+1) -> Cliff(n) because Cliff(n) acts on Vector(n) in an important way.
Second,
The code you wrote can be parametrized on R without any tricks at all, see below. The only troublesome case is 'instance Clifford R'....
Yes, that's exactly where I had trouble. I even played with fundeps but didn't quite get the solution you have.
I'll probably fold your latest addition into the code I added last weekend. I wrote stuff to deal with isomorphisms between algebraic structures eg. Quaternion (Quaternion a) is, in some sense, isomorphic to Matrix (Matrix a), although the isomorphism is non-trivial. Turns out these isomorphisms have a nice application and I'll post something about it eventually.
augustss,
A followup to saying:
In a suitable programming language I'd replace e with something of type e :: Vector(n) -> Cliff(n).
I could have written the code like this from day 1. But I'm interested in eventually figuring out an efficient implementation of this function that doesn't simply map multiples of each basis element in turn and then sum the results. (I think that's what the FFT paper I linked to does but I'm deliberately trying not to read that paper yet...)
The code you wrote can be parameterized on R without any tricks at all, see below. The only troublesome case is 'instance Clifford R'.
It doesn't have to be. Instead of trying to define a general instance, Clifford r r, just write separate instances for any base field you want.
instance Clifford Float Float where
{ i a = a
; e _ _ = undefined
}
Incidentally, what about defining e to return a list of basis vectors? (Assuming I have the terminology right.)
class (Num a) => Clifford r a where
{ i :: r -> a
; e :: [a]
}
instance Clifford Double Double where
{ i a = a
; e = []
}
instance (Clifford r a) => Clifford r (Quaternion a) where
{ i a = Q (i a) 0 0 0
; e = [Q 0 1 0 0, Q 0 0 1 0] ++ [Q 0 0 0 a | a <- e]
}
sigfpe,
So I can see that e :: Vector n -> Cliff n would make sense. But just looking at the code, it looks like 'i' is begging to be 'e 0'. :)
Often when you find symmetries like that it's for some reason.
Dave,
That looks rather neat. :) And it's the right thing to do. The I type was a hack. But does it extend to the MkCliff class in a natural way?
So what this whole thing really needs is dependent types. But not having that, I think we could fake something with type level naturals for getting the type of e right. It really needs to have size n.
That looks rather neat. :) And it's the right thing to do. The I type was a hack. But does it extend to the MkCliff class in a natural way?
Probably, although I haven't tried.
I don't know enough about geometric algebra to judge usage cases, but it seems to me that it would be easier to generate types for Clifford algebras using Template Haskell.
Something along these lines::
cliff 1 r = [t| Complex $r |]
cliff 2 r = [t| Quaternion $r |]
cliff n r = [t| Quaternion $(cliff' (n-2) r) |]
cliff' 1 r = [t| Split $r |]
cliff' 2 r = [t| Matrix $r |]
cliff' n r = [t| Matrix $(cliff (n-2) r) |]
(Unfortunately, GHC won't accept this. Apparently, you can't splice in a type declaration, so you'd have to work directly with the expression tree.)
Going further, it should be possible to write a metaprogram that takes a description of a Clifford algebra and codes up an appropriate optimized implementation.
Dave,
Using Template Haskell isn't always an option. It turns the creation of the type into something static, which isn't always what you want (well, I know next to nothing of Geometric Algebra, so I'm guessing).
Also, using TH gives up static type checking in the same sense as using templates in C++. I don't like that.
Often when you find symmetries like that it's for some reason.
The Clifford algebras are all about symmetries. But overall "e 0" upsets more symmetries than it appears to reveal.
Before I forget again, let me thank sigfpe for this post. I'd heard of geometric algebra, but I'd never been interested enough to find out what it is. Also, it's cool to know what Complex (Complex Float) means.
augustss: Using Template Haskell isn't always an option. It turns the creation of the type into something static, which isn't always what you want [...]
True enough. My suspicion is that most programs using Clifford algebras will know which algebra they want at compile time, but it is possible to defer the choice to run time.
Ignoring the parameterization over the base field, we can write:
cliff :: Int -> (forall a. (Clifford a, Num a) => a -> w) -> w
cliff 1 k = k (undefined :: Cliff1)
cliff 2 k = k (undefined :: Cliff2)
cliff n k = cliff' (n-2) (\(_::a) -> k (undefined :: Quaternion a))
cliff' :: Int -> (forall a. (Clifford a, Num a) => a -> w) -> w
cliff' 1 k = k (undefined :: Cliff1')
cliff' 2 k = k (undefined :: Cliff2')
cliff' n k = cliff (n-2) (\(_::a) -> k (undefined :: Matrix a))
Doing this while parameterizing over the base field is trickier, because we need to construct the dictionaries, but I'm confident it can be done.
The link you gave to my preprint at UNSW does not work.
The reference is:
Paul Leopardi, "A generalized FFT for Clifford algebras", Bulletin of the Belgian Mathematical Society - Simon Stevin, Volume 11, Number 5, 2005, pp. 663-688.
Preprint: UNSW Applied Mathematics Report AMR04/17, March 2004.
See also the GluCat C++ template library.
Hmmm. Very interesting.
I'm coming here from the page about Synthetic Differential Geometry.
This gives at least some partial indications about the relationship of that stuff to Clifford Algebras, at least by implication ... .
Supersymmetry, etc??
Hmmm. Mentions Bott Periodicity. So we are into K-Theory. Maybe get to the Atiyah-Singer Index Theorem.
And we are implementing things in Haskell. Very interesting.
What IS this blog??
@john rood It is what it is :-)
I will want to absorb this. For me, this is a "random find". I was looking for links between dual numbers and synthetic differential geometry on the net and this blog came up.
I will give three hints about myself. I am old--age 62. I just finished Coursera's Exploring Quantum Physics MOOC of Galitski and Clark. I am a coauthor of a paper with Peter Norvig from very long ago (cf Norvig.com).
@john rood I work at the same company as Peter Norvig.
I'm actually not really sure what the relationship is between dual numbers and SDG. Dual numbers have nilpotents because they're an ordinary algebraic structure in a classical framework. SDG has nilpotents because it's in a non-classical logic where for the relevant notion of real number, d^2=0 doesn't imply d=0.
Yes, i figured out that you work at Google. I am a first level link to Peter. (He is very busy and has many of these.) He is my friend actually.
Post a Comment
<< Home