The free vector space on a type, part 1

As I mentioned last time, I want to spend the next few posts talking about quantum algebra. Well, we've got to start somewhere, so let's start with vector spaces.

You probably know what a vector space is. It's what is sounds like: a space of vectors, that you can add together, or multiply by scalars (that is, real numbers, or more generally, elements of some field k). Here's the official definition:
An additive (or Abelian) group is a set with a binary operation called addition, such that
- addition is associative: x+(y+z) = (x+y)+z
- addition is commutative: x+y = y+x
- there is an additive identity: x+0 = x = 0+x
- there are additive inverses: x+(-x) = 0 = (-x)+x
A vector space over a field k is an additive group V, together with an operation k * V -> V called scalar multiplication, such that:
- scalar multiplication distributes over vector addition: a(x+y) = ax+ay
- scalar multiplication distributes over scalar addition: (a+b)x = ax+bx
- associativity: (ab)x = a(bx)
- unit: 1a = a

There are some obvious examples:
- R^2 is a 2-dimensional vector space over the reals R, R^3 is a 3-dimensional vector space. (I'm not going to define dimension quite yet, but hopefully it's intuitively obvious what it means.)
- R^n is an R-vector space for any n.
- Indeed, k^n is a k-vector space for any field k.

Some slightly more interesting examples:
- C is a 2-dimensional vector space over R. (The reason it's more interesting is that of course C possesses additional algebraic structure, beyond the vector space structure.)
- 2*2 matrices over k form a 4-dimensional k-vector space.
- Polynomials in X with coefficients in k form an (infinite dimensional) k-vector space.


If we wanted to code the above definition into Haskell, probably the first idea that would come to mind would be to use type classes:

class AddGrp a where
    add :: a -> a -> a
    zero :: a
    neg :: a -> a -- additive inverse

class (Field k, AddGrp v) => VecSp k v where
    smult :: k -> v -> v -- scalar multiplication

(Type classes similar to these are defined are in the vector-space package.)


For most vector spaces that one encounters in "real life", there is some set of elements, usually obvious, which form a "basis" for the vector space, meaning that all elements can be expressed as linear combinations of basis elements. For example, in R^3, the obvious basis is {(0,0,1), (0,1,0), (1,0,0)}. Any element (x,y,z) of R^3 can be expressed as the linear combination x(1,0,0)+y(0,1,0)+z(0,0,1).

(Mathematicians would want to stress that there are other bases for R^3 that would serve equally well, and indeed, that a significant part of the theory of vector spaces can go through without even talking about bases. However, for our purposes - we want to write code to calculate in vector spaces - then working with a basis is natural.)


Okay, so we want a way to build a vector space from a basis. (More specifically, a k-vector space, for some given field k.) What sorts of things shall we allow as our basis? Well, why not just allow any type, whatsoever:

module Math.Algebras.VectorSpace where

import qualified Data.List as L

data Vect k b = V [(b,k)] deriving (Eq,Ord)

This says that a k-vector space over basis b consists of a linear combination of elements of b. (So the [(b,k)] is to be thought of as a sum, with each (b,k) pair representing a basis element in b multiplied by a scalar coefficient in k.)

For example, we can define the "boring basis" type, which just consists of numbered basis elements:

newtype EBasis = E Int deriving (Eq,Ord)

instance Show EBasis where show (E i) = "e" ++ show i

e i = return (E i) -- don't worry about what "return" is doing here for the moment
e1 = e 1
e2 = e 2
e3 = e 3

Then a typical element of Vect Double EBasis is:

> :load Math.Algebras.VectorSpace
> V [(E 1, 0.5), (E 3, 0.7)]
0.5e1+0.7e3

So of course, the Show instances for EBasis (see above), and Vect k b (not shown) are coming into play here.

How do we know that this is a vector space? Well actually, it's not yet, because we haven't defined the addition and scalar multiplication operations on it. So, without further ado:

infixr 7 *>
infixl 6 <+>

-- |The zero vector
zero :: Vect k b
zero = V []

-- |Addition of vectors
add :: (Ord b, Num k) => Vect k b -> Vect k b -> Vect k b
add (V ts) (V us) = V $ addmerge ts us

-- |Addition of vectors (same as add)
(<+>) :: (Ord b, Num k) => Vect k b -> Vect k b -> Vect k b
(<+>) = add

addmerge ((a,x):ts) ((b,y):us) =
    case compare a b of
    LT -> (a,x) : addmerge ts ((b,y):us)
    EQ -> if x+y == 0 then addmerge ts us else (a,x+y) : addmerge ts us
    GT -> (b,y) : addmerge ((a,x):ts) us
addmerge ts [] = ts
addmerge [] us = us

-- |Negation of vector
neg :: (Num k) => Vect k b -> Vect k b
neg (V ts) = V $ map (\(b,x) -> (b,-x)) ts

-- |Scalar multiplication (on the left)
smultL :: (Num k) => k -> Vect k b -> Vect k b
smultL 0 _ = zero -- V []
smultL k (V ts) = V [(ei,k*xi) | (ei,xi) <- ts]

-- |Same as smultL. Mnemonic is "multiply through (from the left)"
(*>) :: (Num k) => k -> Vect k b -> Vect k b
(*>) = smultL

A few things to mention:
- First, note that we required a Num instance for k. Strictly speaking, as we stated that k is a field, then we should have required a Fractional instance. However, on occasion we are going to break the rules slightly.
- Second, note that for addition, we required an Ord instance for b. We could have defined addition using (++) to concatenate linear combinations - however, the problem with that is that it wouldn't then easily follow that e1+e3 = e3+e1, or that e1+e1 = 2e1. By requiring an Ord instance, we can guarantee that there is a unique normal form in which to express any vector - namely, list the basis elements in order, combine duplicates, remove zero coefficients.
- Finally, note that I didn't define Vect k b as an instance of a vector space type class. That's just because I didn't yet see a reason to.

It will turn out to be useful to have a function that can take an arbitrary element of Vect k b and return a vector in normal form:

-- |Convert an element of Vect k b into normal form. Normal form consists in having the basis elements in ascending order,
-- with no duplicates, and all coefficients non-zero
nf :: (Ord b, Num k) => Vect k b -> Vect k b
nf (V ts) = V $ nf' $ L.sortBy compareFst ts where
    nf' ((b1,x1):(b2,x2):ts) =
        case compare b1 b2 of
        LT -> if x1 == 0 then nf' ((b2,x2):ts) else (b1,x1) : nf' ((b2,x2):ts)
        EQ -> if x1+x2 == 0 then nf' ts else nf' ((b1,x1+x2):ts)
        GT -> error "nf': not pre-sorted"
    nf' [(b,x)] = if x == 0 then [] else [(b,x)]
    nf' [] = []
    compareFst (b1,x1) (b2,x2) = compare b1 b2

Okay, so we ought to check that the Vect k b is a vector space. Let's write some QuickCheck properties:

prop_AddGrp (x,y,z) =
    x <+> (y <+> z) == (x <+> y) <+> z && -- associativity
    x <+> y == y <+> x                 && -- commutativity
    x <+> zero == x                    && -- identity
    x <+> neg x == zero                   -- inverse

prop_VecSp (a,b,x,y,z) =
    prop_AddGrp (x,y,z) &&
    a *> (x <+> y) == a *> x <+> a *> y && -- distributivity through vectors
    (a+b) *> x == a *> x <+> b *> x     && -- distributivity through scalars
    (a*b) *> x == a *> (b *> x)         && -- associativity
    1 *> x == x                            -- unit

instance Arbitrary EBasis where
    arbitrary = do n <- arbitrary :: Gen Int
                   return (E n)

instance Arbitrary Q where
    arbitrary = do n <- arbitrary :: Gen Integer
                   d <- arbitrary :: Gen Integer
                   return (if d == 0 then fromInteger n else fromInteger n / fromInteger d)

instance Arbitrary (Vect Q EBasis) where
    arbitrary = do ts <- arbitrary :: Gen [(EBasis, Q)]
                   return $ nf $ V ts

prop_VecSpQn (a,b,x,y,z) = prop_VecSp (a,b,x,y,z)
    where types = (a,b,x,y,z) :: (Q, Q, Vect Q EBasis, Vect Q EBasis, Vect Q EBasis)

> quickCheck prop_VecSpQn
+++ OK, passed 100 tests.

(I'm using Q instead of R as my field in order to avoid false negatives caused by the fact that arithmetic in Double is not exact.)

So it looks like Vect k b is indeed a vector space. In category theory, it is called the free k-vector space over b. "Free" here means that there are no relations among the basis elements: it will never turn out, for example, that e1 = e2+e3.


Vector spaces of course form a category, specifically an algebraic category (there are other types, as we'll see in due course). The objects in the category are the vector spaces. The arrows or morphisms in the category are the functions between vector spaces which "commute" with the algebraic structure. Specifically, they are the functions f such that:
- f(x+y) = f(x)+f(y)
- f(0) = 0
- f(-x) = -f(x)
- f(a.x) = a.f(x)

Such a function is called linear, the idea being that it preserves lines. This is because it follows from the the conditions that f(a.x+b.y) = a.f(x)+b.f(y) .

We can write a QuickCheck property to check whether a given function is linear:

prop_Linear f (a,x,y) =
    f (x <+> y) == f x <+> f y &&
    f zero == zero             &&
    f (neg x) == neg (f x)     &&
    f (a *> x) == a *> f x

prop_LinearQn f (a,x,y) = prop_Linear f (a,x,y)
    where types = (a,x,y) :: (Q, Vect Q EBasis, Vect Q EBasis)

For example:

> quickCheck (prop_LinearQn (2 *>))
+++ OK, passed 100 tests.

We won't need to use this quite yet, but it's handy to have around.


Now, in category theory we have the concept of a functor, which is a map from one category to another, which commutes with the category structure. Specifically, a functor F consists of a map from the objects of one category to the objects of the other, and from the arrows of one category to the arrows of the other, satisfying:
- F(id_A) = id_F(A)
- F(f . g) = F(f) . F(g) (where dot denotes function composition)

How does this relate to the Functor type class in Haskell? Well, the Haskell type class enables us to declare that a /type constructor/ is a functor. For example, (Vect k) is a type constructor, which acts on a type b to construct another type Vect k b. (Vect k) is indeed a functor, witness the following declaration:

instance Functor (Vect k) where
    fmap f (V ts) = V [(f b, x) | (b,x) <- ts]

This says that if we have a function f on our basis elements, then we can lift it to a function on linear combinations of basis elements in the obvious way.

In mathematics, we would think of the free vector space construction as a functor from Set (the category of sets) to k-Vect (the category of k-vector spaces). In Haskell, we need to think of the (Vect k) construction slightly differently. It operates on types, rather than sets, so the source category is Hask, the category of Haskell types.

What is the relationship between Hask and Set? Well, if we identify a type with the set of values which inhabit it, then we can regard Hask as a subcategory of Set, consisting of those sets and functions which can be represented in Haskell. (That would imply for example that we are restricted to computable functions.)

So (Vect k) is a functor from Hask to the subcategory of k-Vect consisting of vector spaces over sets/types in Hask.

So just to spell it out, the (Vect k) functor:
- Takes an object b in Hask/Set - ie a type, or its set of inhabitants - to an object Vect k b in k-Vect
- Takes an arrow f in Hask (ie a function f :: a -> b), to an arrow (fmap f) :: Vect k a -> Vect k b in k-Vect.

Now, there's just one small fly in the ointment. In order to get equality of vectors to work out right, we wanted to insist that they were expressed in normal form, which meant we needed an Ord instance for b. However, in the Functor instance for (Vect k), Haskell doesn't let us express this constraint, and our fmap is unable to use the Ord instance for b. What this means is that fmap f might return a vector which is not in normal form - so we need to remember to call nf afterwards. For example:

newtype FBasis = F Int deriving (Eq,Ord)

instance Show FBasis where show (F i) = "f" ++ show i

> let f = \(E i) -> F (10 - div i 2)
> let f' = fmap f :: Vect Q EBasis -> Vect Q FBasis
> f' (e1 <+> 2 *> e2 <+> e3)
f10+2f9+f9
> let f'' = nf . fmap f :: Vect Q EBasis -> Vect Q FBasis
> f'' (e1 <+> 2 *> e2 <+> e3)
3f9+f10

So it might be fairer to say that it is the combination of nf and fmap that forms the functor on arrows.

The definition of a functor requires that the target arrow is an arrow in the target category. In this case, the requirement is that it is a linear function, rather than just any function between vector spaces. So let's just check:

> quickCheck (prop_LinearQn f'')
+++ OK, passed 100 tests.

That's enough for now - more next time.