Oleksandr Manzyuk's Blog

Musings of a mathematician and aspiring programmer

Month: October, 2012

Gröbner bases in Haskell: Part II

In the previous post we discussed the representation of variables, monomials, and monomial orderings in a toy EDSL for polynomial computations groebner. In this post, we are going to consider how polynomials are represented in the EDSL and implement the Gröbner basis algorithm.

Polynomials

A term is a monomial times a coefficient (an element of some ground ring r):

data Term r v o = T r (Monomial v o) deriving (Eq, Show)

Terms as well as monomials form a monoid:

instance (Num r, Ord v) => Monoid (Term r v o) where
    mempty = T 1 mempty
    T a m `mappend` T b n = T (a * b) (m `mappend` n)

Polynomials are represented as lists of non-zero terms, ordered in descending order by the their monomials with respect to some monomial order. This makes equality test and extraction of the leading monomial very simple and cheap:

newtype Polynomial r v o = P [Term r v o] deriving Eq

lm :: Polynomial r v o -> Monomial v o
lm (P ((T _ m):_)) = m
lm (P [])          = error "lm: zero polynomial"

The relative order of terms does not change if the polynomial is multiplied by a term:

(*^) :: (Num r, Ord v) => Term r v o -> Polynomial r v o -> Polynomial r v o
u *^ P vs = P [ u `mappend` v | v <- vs ]

Polynomials have degree: the degree of the zero polynomial is usually defined to be -1, and the degree of a non-zero polynomial is the maximum of the degrees of its monomials:

instance Ord v => HasDegree (Polynomial r v o) where
    degree (P []) = -1
    degree (P ts) = maximum [ degree m | T _ m <- ts ]

We are trying to make the display of polynomials as close to the mathematical notation as possible. Because we don’t know what the ground ring r can be, we apply some heuristics:

instance (Eq r, Show r, Num r, Ord v, Show v) => Show (Polynomial r v o) where
    show (P [])     = "0"
    show (P (t:ts)) = showHead t ++ showTail ts
        where
          showHead (T c m) = prefix ++ show m
              where
                prefix = case c of
                           1  -> ""
                           -1 -> "-"
                           _  -> show c
          showTerm (T c m) = prefix ++ show m
              where
                prefix = case signum c of
                           1  -> '+':a
                           -1 -> '-':a
                           _  -> "(" ++ show c ++ ")"
                a = if abs c == 1 then "" else show (abs c)
          showTail = concatMap showTerm

Arithmetic operations on polynomials are defined to preserve the invariant of the representation of polynomials:

instance (Eq r, Num r, Ord v, Show v, Ord (Monomial v o))
    => Num (Polynomial r v o) where
    f@(P (u@(T a m):us)) + g@(P (v@(T b n):vs))
        | m == n && a + b /= 0
        = let P ws = P us + P vs in P $ T (a + b) m:ws
        | m == n && a + b == 0
        = P us + P vs
        | m < n
        = let P ws = f + P vs in P $ v:ws
        | otherwise
        = let P ws = P us + g in P $ u:ws
    f + P [] = f
    P [] + g = g

    P (u:us) * P (v:vs)
        = let P ws = P us * P vs + u *^ P vs + v *^ P us
          in P $ (u `mappend` v):ws
    _ * P [] = P []
    P [] * _ = P []

    negate (P ts) = P $ [ T (negate a) m | T a m <- ts ]
    -- Inclusion of 'abs' and 'signum' into 'Num' was a stupid idea.
    abs _ = error "abs is undefined for polynomials"
    signum _ = error "signum is undefined for polynomials"
    fromInteger = constant . fromInteger

Constants can be viewed as polynomials (of degree 0 unless the constant is 0):

constant :: (Eq r, Num r, Ord v) => r -> Polynomial r v o
constant 0 = P []
constant c = P [T c mempty]

Variables can be viewed as polynomials of degree 1:

variable :: (Num r, Eq v) => v -> Polynomial r v o
variable x = P [T 1 (inject x)]

Suppose that f and g are two polynomials over a field r, and let x^\alpha and x^\beta be the leading monomials of f and g respectively. Let x^\gamma be the least common multiple of x^\alpha and x^\beta. Then the s-polynomial of f and g is defined to be

\displaystyle \mathrm{spoly}(f, g) = x^{\gamma-\alpha}f - \frac{\mathrm{lc}(f)}{\mathrm{lc}(g)}x^{\gamma-\beta}g

where \mathrm{lc} denotes the leading coefficient of a polynomial. In Haskell:

spoly :: (Eq r, Fractional r, Ord v, Show v, Ord (Monomial v o))
      => Polynomial r v o -> Polynomial r v o -> Polynomial r v o
spoly f@(P (u@(T a m):us)) g@(P (v@(T b n):vs)) = n' *^ f - m' *^ g
    where
      n' = T 1       (complement m n)
      m' = T (a / b) (complement n m)

Normal Forms and Gröbner Bases

A normal form of a polynomial f with respect to a list s of polynomials is essentially the remainder from multivariate division of f by polynomials from the list s. It is computed using Buchberger’s algorithm. Instead of explaining the algorithm in words, I’ll let the code speak for itself. The following definition is in fact very close to the pseudo-code that can be found, for example, in A Singular Introduction to Commutative Algebra:

nf :: (Eq r, Fractional r, Ord v, Show v, Ord (Monomial v o))
   => Polynomial r v o -> [Polynomial r v o] -> Polynomial r v o
nf f s = go f
    where
      go h | h == 0      = 0
           | []    <- s' = h
           | (g:_) <- s' = go (spoly h g)
           where
             s' = [g | g <- s, lm h `isDivisibleBy` lm g]

The function groebner implements the Gröbner basis algorithm. It takes a list of generators of an ideal and returns a Gröbner basis of that ideal:

groebner :: (Eq r, Fractional r, Ord v, Show v, Ord (Monomial v o))
         => [Polynomial r v o] -> [Polynomial r v o]
groebner i = go i ps
    where
      ps = [(f, g) | f <- i, g <- i, f /= g]
      go s [] = s
      go s ps@((f, g):ps')
          | h == 0    = go s ps'
          | otherwise = go (h:s) (ps' ++ [(h, f) | f <- s])
          where
            h = nf (spoly f g) s

Product Criterion

The product criterion allows us to decrease the number of pairs that have to be considered by the Groebner basis algorithm. The criterion says that if the least common multiple of the leading monomials of f and g is their product, then the s-polynomial of f and g reduces to 0 with respect to the set \{f, g\}, and hence the pair (f, g) can be dropped. We implement a function pc that tests if two polynomials f and g satisfy the product criterion as follows:

pc f g = null (variables (lm f) `intersect` variables (lm g))

Note that instead of computing the least common multiple of the leading monomials of f and g we check if the sets of variables occurring in each of the monomials are disjoint.

The optimized version of the Gröbner basis algorithm reads as follows:

groebner :: (Eq r, Fractional r, Ord v, Show v, Ord (Monomial v o))
         => [Polynomial r v o] -> [Polynomial r v o]
groebner i = go i ps
    where
      ps = [(f, g) | f <- i, g <- i, f /= g, not (pc f g)]
      go s [] = s
      go s ps@((f, g):ps')
          | h == 0    = go s ps'
          | otherwise = go (h:s) (ps' ++ [(h, f) | f <- s, not (pc h f)])
          where
            h = nf (spoly f g) s

Declaring Variables

Having to define an enumeration type and to write an Enumerable instance for it in order to declare variables is tedious. Here we address this problem.

First, let us define the sum of types a and b:

data a :<: b = Inl a | Inr b deriving (Eq, Ord)
infixr 6 :<:

instance (Show a, Show b) => Show (a :<: b) where
    show (Inl x) = show x
    show (Inr x) = show x

It will become clear shortly why we have chosen to denote the sum by :<:. If both a and b are enumerable, then so is the sum of a and b:

instance (Enumerable a, Enumerable b) => Enumerable (a :<: b) where
    enumerate = map Inl enumerate ++ map Inr enumerate

Note that the elements of a are enumerated before the elements of b.

We can now define each variable as a singleton type:

data X = X
data Y = Y
data Z = Z

and join these types into X :<: Y :<: Z. Each of the types X, Y, Z is Enumerable:

instance Enumerable X where enumerate = [X]
instance Enumerable Y where enumerate = [Y]
instance Enumerable Z where enumerate = [Z]

Hence, the type X :<: Y :<: Z is Enumerable too. It is isomorphic to

data XYZ = X | Y | Z

However, the elements of X :<: Y :<: Z are somewhat unwieldy to write: Inl X, Inr (Inl Y), and Inr (Inr Z). We solve this problem by the trick used in Data types à la carte. We introduce the following typeclass that expresses the fact that a type a is a subtype of a type b:

class Sub a b where
    inj :: a -> b

Instead of writing the injections using Inl and Inr, the injections will be inferred using this typeclass. The Sub typeclass has only three instances:

a is a subtype of a:

instance Sub a a where
    inj = id

a is also a subtype of the sum a :<: b:

instance Sub a (a :<: b) where
    inj = Inl

Finally, if a is a subtype of c, then a is also a subtype of the sum b :<: c:

instance Sub a c => Sub a (b :<: c) where
    inj = Inr . inj

These instances require quite a few GHC extensions including OverlappingInstances (in fact, we have implicitly used many GHC extensions above), but as with data types à la carte, this shouldn’t result in an unexpected behaviour, provided that user never explicitly nests sums.

The following function allows us to view a variable from a set v as a polynomial over any wider set of variables w:

var :: (Sub v w, Ord (Monomial w o), Num r, Eq w) => v -> Polynomial r w o
var = variable . inj

We can now introduce

x = var X
y = var Y
    z = var Z

Thus, x is a polynomial over any ring r, over any set of variables w containing X, with respect to any monomial ordering:

*Main> :t x
x :: (Eq w, Num r, Ord (Monomial w o), Sub X w) => Polynomial r w o

Ditto for y and z.

Consequently, for example, x * y + z is a polynomial over any ring r, over any set of variables w containing X, Y, and Z, in particular it is a polynomial over the set of variables X :<: Y :<: Z.

Writing separate definitions for X, Y, Z and Enumerable instances is still tedious. Fortunately, we can write a Template Haskell macro allowing us to write

$(defineVariables ["X", "Y", "Z"])

See Variable.hs for details.

Example

Consider the ideal generated by polynomials x^{10} + x^9y^2 and y^8 - x^2y^7:

ideal :: Ord (Monomial (X :<: Y) o) => [Polynomial Rational (X :<: Y) o]
ideal = [x ^ 10 + x ^ 9 * y ^ 2, y ^ 8 - x ^ 2 * y ^ 7]

Note that ideal is polymorphic in the monomial ordering. Let basis be the Gröbner basis of ideal:

basis :: Ord (Monomial (X :<: Y) o) => [Polynomial Rational (X :<: Y) o]
basis = groebner ideal

We can now compute basis with respect to different monomial orderings:

*Main> basis :: [Polynomial Rational (X :<: Y) Lex]
[Y^15-Y^12,XY^12+Y^14,XY^13+Y^12,X^10+X^9Y^2,-X^2Y^7+Y^8]
*Main> basis :: [Polynomial Rational (X :<: Y) RevLex]
[X^16+X^13,-X^30-X^27,-X^13Y+X^15,-X^14Y-X^13,X^9Y^2+X^10,Y^8-X^2Y^7]
*Main> basis :: [Polynomial Rational (X :<: Y) DegLex]
[Y^14+XY^12,Y^18-X^4Y^13,XY^13+X^2Y^11,XY^17-X^11Y^6,-X^13+XY^12,X^12Y+X^3Y^10,X^9Y^2+X^10,-X^2Y^7+Y^8]
*Main> basis :: [Polynomial Rational (X :<: Y) DegRevLex]
[Y^14+XY^12,Y^18-X^4Y^13,XY^13+X^2Y^11,XY^17-X^11Y^6,-X^13+XY^12,X^12Y+X^3Y^10,X^9Y^2+X^10,-X^2Y^7+Y^8]

Advertisements

Gröbner bases in Haskell: Part I

Introduction

Although I am primarily a category theorist (my PhD research was at the border between category theory and homological algebra), I also studied algebraic geometry and computer algebra for my master degree at the University of Kaiserslautern. Computer algebra is a diverse subject encompassing many fields of symbolic computations. Our lectures in Kaiserslautern focused on computations in polynomial rings with an eye towards applications to algebraic geometry.

Functional programming languages are a particularly good fit for symbolic computations. My goal in this blog post is to illustrate this by showing you an implementation in Haskell of Buchberger’s algorithm for computing Gröbner bases of polynomial ideals.

Haskell is also a great host language for embedded domain-specific languages. To illustrate this, we are going to implement Buchberger’s algorithm as part of a tiny EDSL for polynomial computations.

Before we dive in, I would like to make a couple of remarks about the desiderata that guided the implementation.

First, my goal was to design a toy EDSL, not a fully-fledged computer algebra system, with the syntax that is as close to the mathematical notation as possible, possibly at the cost of efficiency. Don’t expect this implementation to be competitive with the advanced computer algebra systems like Singular, Macaulay, Maple etc.

Second, the whole point of using Haskell, and not say Common Lisp, as the host language is to leverage Haskell’s highly expressive type system. For example, in the computer algebra system Singular, developed in Kaiserslautern, polynomial rings are described by three pieces of information: the ground field, list of variables, and monomial ordering. In my EDSL this information is encoded in types, so that polynomials over different fields, or in different sets of variables, or of different monomial orderings are values of different types and can’t be accidentally mixed in programs. This particular desideratum has influenced many design choices we are about to discuss.

Today we are going to look at how variables, monomials, and monomial orderings are represented, leaving the representation of polynomials and actual implementation of the Gröbner basis algorithm for the next post. The complete code of the EDSL is available at GitHub.

Variables

Somewhat unconventionally, I have chosen to represent variables by enumeration types. For example, the set of three variables X, Y, Z is represented by the following data type:

data XYZ = X | Y | Z deriving (Eq, Ord, Show)

Every constructor is nullary and simply names a variable. If you want to compute with polynomials in n variables, you first need to define an enumeration type with n constructors, e.g., if n = 5:

data Vars = X1 | X2 | X3 | X4 | X5 deriving (Eq, Ord, Show)

Defining enumeration types by hand is arguably tedious and ugly (having to invent names for enumeration types is particularly obnoxious). We shall see later how we can make declaring variables more pleasant using some Template Haskell.

Monomials

A monomial over a finite set of variables represented by an enumeration type v is represented internally by a Map from variables to integers (exponents):

newtype Monomial v o = M (Map v Int) deriving Eq

The type Monomial is a phantom type: the type variable o does not appear on the right hand side of the definition of Monomial. We are going to use o later as a tag allowing us to define different instances of the Ord typeclass (i.e., different monomial orderings) on Monomial.

We have a number of helper functions to construct monomials. A variable can be viewed as a monomial:

inject :: Eq v => v -> Monomial v o
inject x = M $ Map.singleton x 1

We can convert a monomial to a list of variable-exponent pairs, and we can build a monomial from such a list:

toList :: Ord v => Monomial v o -> [(v, Int)]
toList (M m) = [ p | p@(x, n) <- Map.toList m, n /= 0 ]

fromList :: Ord v => [(v, Int)] -> Monomial v o
fromList xs = M $ Map.fromList [ p | p@(x, n) <- xs, n /= 0 ]

Note that the variables with zero exponents are dropped. We can look up the exponent of a given variable in a monomial:

exponent :: Ord v => v -> Monomial v o -> Int
exponent x (M m) = fromMaybe 0 (Map.lookup x m)

We can also collect the variables occurring in a monomial with non-zero exponents:

variables :: Ord v => Monomial v o -> [v]
variables = map fst . toList

Monomials are shown as power products of variables:

instance (Ord v, Show v) => Show (Monomial v o) where
    show m
        | null support
        = "1"
        | otherwise
        = concat [ show x ++ suffix
                 | (x, n) <- support
                 , let suffix = if n == 1
                                then ""
                                else "^" ++ show n
                 ]
        where
          support = toList m

Monomials over a set of variables v are naturally a monoid:

instance Ord v => Monoid (Monomial v o) where
    mempty = M Map.empty
    M a `mappend` M b = M $ Map.unionWith (+) a b

Monomials have degree. Because polynomials also have degree, it is convenient to have an overloaded function degree:

class HasDegree a where
  degree :: a -> Int

The degree of a monomial is the sum of the exponents of its variables:

instance Ord v => HasDegree (Monomial v o) where
    degree (M m) = Map.fold (+) 0 m

We can test whether one monomial is divisible by another:

isDivisibleBy :: Ord v => Monomial v o -> Monomial v o -> Bool
isDivisibleBy (M a) (M b) = Map.isSubmapOfBy (<=) b a

We can divide one monomial by another:

div :: Ord v => Monomial v o -> Monomial v o -> Monomial v o
div (M a) (M b) = M $ Map.differenceWith sub a b
    where
      sub x y | x > y     = Just (x - y)
              | otherwise = Nothing

We are going to use the function div only when one monomial is known to be divisible by the other, but the above definition gives a plausible answer also when this is not the case.

The least common multiple of monomials is defined as follows:

lcm :: Ord v => Monomial v o -> Monomial v o -> Monomial v o
lcm (M a) (M b) = M $ Map.unionWith max a b

Finally, complement m n computes the product of factors in n that are missing in m:

complement :: Ord v => Monomial v o -> Monomial v o -> Monomial v o
complement m n = lcm m n `div` m

Monomial Orderings

A monomial ordering on the set of monomials is a total ordering that is compatible with multiplication of monomials, i.e., m_1 \le m_2 implies m_1n \le m_2n.

We equip monomials Monomial v o over the set of variables v with different orderings by supplying different tags o:

instance (Show v, Enumerable v) => Ord (Monomial v Lex) where
    (<=) = lex
instance (Show v, Enumerable v) => Ord (Monomial v RevLex) where
    (<=) = revlex
instance (Show v, Enumerable v) => Ord (Monomial v DegLex) where
    (<=) = deglex
instance (Show v, Enumerable v) => Ord (Monomial v DegRevLex) where
    (<=) = degrevlex

Here Lex, RevLex, DegLex, and DegRevLex are empty data types:

data Lex         -- Lexicographic ordering
data RevLex      -- Reverse lexicographic ordering
data DegLex      -- Degree lexicographic ordering
data DegRevLex   -- Reverse degree lexicographic ordering

that are used as tags so that we can define different Ord instances on the Monomial type. Instead of making Monomial a phantom type, we could also define newtype wrappers around Map v Int, but then the definition of polynomials would have to become more involved.

The definitions of orderings are written in a slightly unintuitive style because they define (<=), not (>) or (<) as is customary. This is necessary because a minimal instance declaration of Ord requires either compare or (<=). In particular, if we define only (>), then the default implementation of (<=) isn’t in terms of (>) but in terms of compare, which in turn by default is defined in terms of (<=), leading to an infinite loop.

lex' :: (Ord v, Show v) => Monomial v o -> Monomial v o -> [v] -> Bool
lex' a b []     = True
lex' a b (x:xs) = exponent x a <= exponent x b
                  && (exponent x a /= exponent x b || lex' a b xs)

lex, revlex, deglex, degrevlex :: (Enumerable v, Show v)
                               => Monomial v o -> Monomial v o -> Bool
lex       a b = lex' a b enumerate
revlex    a b = lex' a b (reverse enumerate)
deglex    a b = degree a <= degree b
                && (degree a /= degree b || a `lex` b)
degrevlex a b = degree a <= degree b
                && (degree a /= degree b || b `revlex` a)

The definitions of orderings rely on the order of variables and the knowledge which variables can occur in monomials. We encode this information in the typeclass Enumerable:

class Ord a => Enumerable a where
    enumerate :: [a]

We have added the Ord constraint to the context in order to save some typing (and because it makes sense). For each type v of variables, we have to define an instance of Enumerable. For example:

data XYZ = X | Y | Z deriving (Eq, Ord, Show)

instance Enumerable XYZ where enumerate = [X, Y, Z]

Summary

That’s it for this time. Thanks for reading this far! In this blog post we have discussed how variables, monomials, and monomial orderings are represented in our EDSL. In the next post we shall consider how polynomials are represented in the EDSL and implement the Gröbner basis algorithm.