Sunday, 18 September 2011

Commutative Algebra and Algebraic Geometry

Last time we saw how to create variables for use in polynomial arithmetic. This time I want to look at some of the things we can do next.

First, let's define the variables we are going to use:

> :l Math.CommutativeAlgebra.GroebnerBasis
> let [t,u,v,x,y,z,x',y',z'] = map glexvar ["t","u","v","x","y","z","x'","y'","z'"]

So now we can do arithmetic in the polynomial ring Q[t,u,v,x,y,z,x',y',z']. For example:

> (x+y)^2

The branch of mathematics dealing with the theory of polynomial rings is called commutative algebra, and it was "invented" mainly in support of algebraic geometry. Algebraic geometry is roughly the study of the curves, surfaces, etc that arise as the solution sets of polynomial equations. For example, the solution-set of the equation x^2+y^2=1 is the unit circle.

If we are given any polynomial equation f = g, then we can rewrite it more conveniently as f-g = 0. In other words, we only need to track individual polynomials, rather than pairs of polynomials. Call the solution set of f = 0 the zero-set of f.

Sometimes we're interested in the intersection of two or more curves, surfaces, etc. For example, it is well known that the hyperbola, parabola and ellipse all arise as "conic sections" - that is, as the intersection of a cone with a plane. So define the zero-set of a collection (or system) of polynomials to be the set of points which are zeros of all the polynomials simultaneously. For example, the zero-set of the system [x^2+y^2-z^2, z-1] is the unit circle x^2+y^2=1 situated on the plane z=1.

Okay, so how can commutative algebra help us to investigate curves and surfaces? Is there a way for us to "do geometry by doing algebra"? Well, first, what does "doing geometry" consist of? Well, at least some of the following:
- Looking at the shapes of curves and surfaces
- Looking at intersections, unions, differences and products of curves and surfaces
- Looking at when curves or surfaces can be mapped into or onto other curves or surfaces
- Looking at when two different curves or surfaces are equivalent, in some sense (for example, topologically equivalent)

(That phrase "curves and surfaces" is not only clumsy but also inaccurate, so from now on I'll use the proper term, "variety", for the zero-set of a system of polynomials, whether it's a set of isolated points, a curve, a surface, some higher dimensional thing, or a combination of some of the preceding.)

So can we do all those things using algebra? Well, let's have a go.

Let's start by looking at intersections and unions of varieties (remember, that's just the fancy name for curves, surfaces, etc.).

Well, we've already seen how to do intersections. If a variety V1 is defined by a system of polynomials [], and a variety V2 is defined by [], then their intersection is defined by the system [,] - the zero-set of both sets of polynomials simultaneously. We'll call this the "sum" of the systems of polynomials. (Note to the cognoscenti: yes, I'm really talking about ideals here.)

sumI fs gs = gb (fs ++ gs)

Don't worry too much about what that "gb" (Groebner basis) call is doing. Let's just say that it's choosing the best way to represent the system of polynomials. For example:

> sumI [x^2+y^2-z^2] [z-1]

Notice how the gb call has caused the first polynomial to be simplified slightly. The same variety might arise as the zero-set of many different systems of polynomials. That's something that we're going to need to look into - but later.

Okay, so what about unions of varieties. So suppose V1 is defined by [], V2 is defined by []. A point in their union is in either V1 or V2, so it is in the zero-set of either [] or []. So how about multiplying the polynomials together in pairs. That is, let's look at the system [fi*gj | fi <- fs, gj <- gs]. Call the zero-set of this system V. Then clearly, any point in either V1 or V2 is in V, since we then know that either all the fs or all the gs vanish at that point, and hence so do all the products. Conversely, suppose that some point is not in the union of V1 and V2. Then there must exist some fi, and some gj, which are non-zero at that point. Hence there is an fi*gj which is non-zero, so the point is not in V.

This construction is called, naturally enough, the product of the systems of polynomials.

productI fs gs = gb [f * g | f <- fs, g <- gs]

> productI [x^2+y^2-z^2] [z-1]

Just in case you're still a little unsure, let's confirm that a few arbitrary points in the union are in the zero-set of this polynomial:

> eval (x^2*z+y^2*z-z^3-x^2-y^2+z^2) [(x,100),(y,-100),(z,1)]
> eval (x^2*z+y^2*z-z^3-x^2-y^2+z^2) [(x,3),(y,4),(z,5)]

The first expression evaluates the polynomial at the point (100,-100,1), an arbitrary point on the plane z=1. The second evaluates at (3,4,5), an arbitrary point on the cone x^2+y^2=z^2. Both points are in the zero-set of our product polynomial.

Since we're in the neighbourhood, let's have a look at the other conic sections. First, let's rotate our coordinate system by 45 degrees, using the substitution x'=x+z, z'=z-x. (Okay, so this also scales - to save us having to handle a sqrt 2 factor.)

> let cone' = subst (x^2+y^2-z^2) [(x,(x'-z')/2),(y,y'),(z,(x'+z')/2)]
> cone'

In these coordinates, the intersection of the cone with the plane z'=1 is the parabola x'=y'^2:

> sumI [cone'] [z'-1]

Alternatively, the intersection with the plane y'=1 is the hyperbola x'z'=1:

> sumI [cone'] [y'-1]

Okay, so we've made a start on seeing how to do geometry by doing algebra, by looking at unions and intersections of varieties. There's still plenty more to do. We mustn't forget that we have some unfinished business: we need to understand when different polynomial systems can define the same variety, and in what sense the gb (Groebner basis) function finds the "best" representation. That will have to wait for another time.

Incidentally, for the eval and subst functions that I used above, you will need to take the new release HaskellForMaths v0.4.0. In this release I also removed the older commutative algebra modules, so I revved the minor version number.

Sunday, 4 September 2011

Commutative Algebra in Haskell, part 1

Once again, it's been a little while since my last post, and once again, my excuse is partly that I've been too busy writing code.

I've just uploaded a new release, HaskellForMaths 0.3.4, which contains the following new modules:

Math.Core.Utils - this is a collection of utility functions used throughout the rest of the library. I've belatedly decided that it's better to put them all in one place rather than scattered here and there throughout other modules.

Math.Core.Field - this provides new, more efficient implementations of several finite fields. There already were implementations of these finite fields, in the Math.Algebra.Field.Base and ...Extension modules, as discussed here and here. However, that code was written to make the maths clear, rather than for speed. This new module is about speed. For the prime power fields in particular (eg F4, F8, F9), these implementations are significantly faster.

Math.Combinatorics.Matroid - Matroids are a kind of combinatorial abstraction of the concept of linear independence. They're something that I heard about years ago - both of my favourite combinatorics books have brief introductions - but I never bothered to follow up. Well anyway, so something finally piqued my curiosity, and I got Oxley's Matroid Theory. It turned out to be really interesting stuff, and this module is pretty much a translation of a large part of that book into Haskell code, written as I taught myself all about matroids.

Math.CommutativeAlgebra.Polynomial - Although I hadn't yet got around to discussing them in the blog, HaskellForMaths has always had modules for working with multivariate polynomials, namely Math.Algebra.Commutative.Monomial and ...MPoly. However, these were some of the earliest code I wrote, before my more recent free vector space and algebra code. So I saw an opportunity to simplify and improve this code, by building it on top of the free vector space code. Also, I'm trying to rationalise the module naming convention in HaskellForMaths, to more closely follow the categories used in or . In the long run, I expect this module to supercede the older modules.

Math.CommutativeAlgebra.GroebnerBasis - Again, there was already code for Groebner bases in Math.Algebra.Commutative.GBasis. This is pretty much the same code, ported to the new polynomial implementation, but I've also begun to build on this, with code to find the sum, product, intersection, and quotient of ideals.

So the matroid code was just new code that I wrote while teaching myself some new maths. But most of the other code comes from an ambition to organise and simplify the HaskellForMaths library. I've also been trying to improve the documentation.

My ultimate ambition is to get more people using the library. To do that, the structure of the library needs to be clearer, the documentation needs to be better, and I need to explain how to use it. So I thought I'd start by explaining how to use the new commutative algebra modules.

(So this is a bit of a digression from the series on quantum algebra that I've been doing the last few months. However, in terms of the cumulative nature of maths, it's probably better to do this first.)

Okay, so suppose we want to do some polynomial arithmetic. Well, first we need to create some variables to work with. How do we do that?

First, decide on a monomial ordering - that is, we need to decide in what order monomials are to be listed within a polynomial. For the moment, let's use "graded lexicographic" or Glex order. This says that you should put monomials of higher degree before those of lower degree (eg y^3 before x^2), and if two monomials have the same degree, you should use lexicographic (dictionary) order (eg xyz before y^3).

Next, decide on a field to work over. Most often, we'll want to work over Q, the rationals.

Then, our variables themselves can be of any Haskell type - but there are usually only two sensible choices:

The easiest way is to use String as the type for our variables.

Then we could make some variables like this:

> :l Math.CommutativeAlgebra.Polynomial
> let [x,y,z] = map glexvar ["x","y","z"]

And then we can do polynomial arithmetic:

> (x+y+z)^3

If we want to use any other field besides Q, then we will have to use a type annotation to tell the compiler which field we're working over:

> let [x,y,z] = map var ["x","y","z"] :: [GlexPoly F3 String]
> (x+y+z)^3

The alternative to using String for our variables is to define our own type. For example

data Var = X | Y | Z | W deriving (Eq,Ord)

instance Show Var where
    show X = "x"
    show Y = "y"
    show Z = "z"
    show W = "w"

[x,y,z,w] = map glexvar [X,Y,Z,W]

So there you have it - now you can do polynomial arithmetic in Haskell.

So how does it work?

Well, fundamentally, k-polynomials are a free k-vector space on the basis of monomials. So we define a type to implement monomials:

data MonImpl v = M Int [(v,Int)] deriving (Eq)
-- The initial Int is the degree of the monomial. Storing it speeds up equality tests and comparisons

instance Show v => Show (MonImpl v) where
    show (M _ []) = "1"
    show (M _ xis) = concatMap (\(x,i) -> if i==1 then showVar x else showVar x ++ "^" ++ show i) xis
        where showVar x = filter ( /= '"' ) (show x) -- in case v == String

Notice that our monomial implementation is polymorphic in v, the type of the variables.

Next, monomials form a monoid, so we make them an instance of Mon (the HaskellForMaths class for monoids):

instance (Ord v) => Mon (MonImpl v) where
    munit = M 0 []
    mmult (M si xis) (M sj yjs) = M (si+sj) $ addmerge xis yjs

In principle, all we need to do now is define an Ord instance, and then an Algebra instance, using the monoid algebra construction.

However, for reasons that will become clear in future postings, we want to be able to work with various different orderings on monomials, such as Lex, Glex, or Grevlex. So we provide various newtype wrappers around this basic monomial implementation. Here's the code for the Glex ordering that we used above:

newtype Glex v = Glex (MonImpl v) deriving (Eq, Mon) -- GeneralizedNewtypeDeriving

instance Show v => Show (Glex v) where
    show (Glex m) = show m

instance Ord v => Ord (Glex v) where
    compare (Glex (M si xis)) (Glex (M sj yjs)) =
        compare (-si, [(x,-i) | (x,i) <- xis]) (-sj, [(y,-j) | (y,j) <- yjs])

type GlexPoly k v = Vect k (Glex v)

glexvar :: v -> GlexPoly Q v
glexvar v = return $ Glex $ M 1 [(v,1)]

instance (Num k, Ord v, Show v) => Algebra k (Glex v) where
    unit x = x *> return munit
    mult xy = nf $ fmap (\(a,b) -> a `mmult` b) xy

We also have similar newtypes for Lex and Grevlex orderings, which I'll discuss another time.

And that's pretty much it. Now that we have an instance of Algebra k (Glex v), we get a Num instance for free, so we get +, -, *, and fromInteger. That means we can enter expressions like the following:

> (2*x^2-y*z)^2

Note that division is not supported: you can't write x/y, for example. However, as a convenience, I have defined a partial instance of Fractional, which does let you divide by scalars. That means that it's okay to write x/2, for example.

Next time, some more things you can do with commutative algebra.