Oleksandr Manzyuk's Blog

Musings of a mathematician and aspiring programmer

Transducers are monoid homomorphisms

Well, that’s not quite right: monoidal transducers, which are closely related to transducers, are in bijection with monoid homomorphisms between free monoids. This wouldn’t make a catchy title, however, would it?

The following is my humble exploration of the ideas from Franklin Chen’s post, in which he is trying to give a more typeful account of transducers, a new abstraction coming to Clojure.

The type of transducers, or reducing function transformers, from a to b is defined as

type Transducer a b = forall r. Reducer a r -> Reducer b r

where

type Reducer a r = r -> a -> r

is a type of reducers, or reducing functions. Notice the universal quantification over r in the definition of the transducer type. What properties of transducers does it imply?

Let us replace the definition of Reducer with an isomorphic one:

type Reducer a r = a -> Endo r

Endo r is a newtype wrapper around r -> r defined in Data.Monoid, which allows us to give r -> r the structure of a monoid:

instance Monoid (Endo r) where
  mempty = Endo id
  Endo f `mappend` Endo g = Endo (f . g)

We could replace Endo r in the definition of the transducer type

type Transducer a b = forall r. (a -> Endo r) -> (b -> Endo r)

with an arbitrary monoid m and consider monoidal transducers:

type MonoidalTransducer a b = forall m. Monoid m => (a -> m) -> (b -> m)

Clearly, a monoidal transducer gives rise to a transducer:

psi :: MonoidalTransducer a b -> Transducer a b
psi h = h

Remarkably, we can go in the opposite direction as well. The trick is to embed a monoid m into Endo m using the Cayley representation:

rep :: Monoid m => m -> Endo m
rep = Endo . mappend

The function rep is a monoid homomorphism. Furthermore, it is split injective with a splitting

abs :: Monoid m => Endo m -> m
abs (Endo f) = f mempty

That is, abs . rep is the identity function. This allows us to define

phi :: Transducer a b -> MonoidalTransducer a b
phi t f = abs . t (rep . f)

I don’t know if phi and psi are mutually inverse. I thought I had a proof that the composite phi . psi is the identity function, but it relied on the free theorem for the type MonoidalTransducer a b, which requires that abs be a monoid homomorphism, which it is not.

Monoidal transducers are easier to understand than transducers. They can be interpreted categorically as follows.

We have two categories:

  • \mathbf{Set}, the category of sets, whose objects are sets and whose morphisms are ordinary functions, and
  • \mathbf{Mon}, the category of monoids, whose objects are monoids and whose morphisms are monoid homomorphisms.

We have two functors between these categories:

  • U : \mathbf{Mon}\to\mathbf{Set}, the forgetful functor, taking a monoid to its underlying set, and
  • F : \mathbf{Set}\to\mathbf{Mon}, the free monoid functor, taking a set X to a free monoid generated by X, which is nothing but the set of lists over X with the standard structure of a monoid.

The functors U and F are adjoint, i.e., there is a natural bijection between the set of functions \mathbf{Set}(X, UY) and the set of monoid homomorphisms \mathbf{Mon}(FX, Y), for each set X and monoid Y.

Monoidal transducers from a set X to a set Y are interpreted as natural transformations between the functors \mathbf{Set}(X, U(-)) and \mathbf{Set}(Y, U(-)) from the category \mathbf{Mon} to \mathbf{Set}. By adjointness, we can replace the functors \mathbf{Set}(X, U(-)) and \mathbf{Set}(Y, U(-)) with isomorphic functors \mathbf{Mon}(FX, -) and \mathbf{Mon}(FY, -). Therefore, the set of monoidal transducers is isomorphic to the set of natural transformations from \mathbf{Mon}(FX, -) to \mathbf{Mon}(FY, -), which by Yoneda Lemma is isomorphic to the set \mathbf{Mon}(FY, FX), i.e., the set of monoid homomorphism from the free monoid FY to the free monoid FX.

If one were so inclined, one could say that the category of monoidal transducers, whose objects are sets and whose morphisms are monoidal transducers, is equivalent (in fact, isomorphic) to the category of free monoids, a full subcategory of the category of monoids. The latter is isomorphic to the Kleisli category of the list monad: every monoid homomorphism FY\to FX is uniquely determined by its restriction to the set of generators Y, i.e., by the function Y\to FX, and composition of monoid homomorphisms translates into Kleisli composition.

To sum up, monoidal transducers from a to b are in bijection with monoid homomorphisms from [b] to [a], which are in turn in bijection with functions b -> [a]. The bijections are given explicitly by

chi :: MonoidalTransducer a b -> b -> [a]
chi f = f return

rho :: (b -> [a]) -> MonoidalTransducer a b
rho f k = mconcat . map k . f

The monoid homomorphism associated with a function f :: b -> [a] is concatMap f :: [b] -> [a]. Note that the functions

mapping :: (a -> b) -> MonoidalTransducer b a
mapping f k = k . f

filtering :: (a -> Bool) -> MonoidalTransducer a a
filtering p k x = if p x then k x else mempty

give rise to the functions map and filter on lists, which are monoid homomorphisms. Conversely, it is impossible to define a function

taking :: Int -> MonoidalTransducer a a

that would give rise to the function take because take is not a monoid homomorphism.

From Object Algebras to Finally Tagless Interpreters

Last Saturday I gave a talk at kiev::fprog meetup. I really enjoyed the event, and taking the opportunity I would like to thank Vladimir Kirillov for organizing it.

The slides of the talk can be found here, but I thought I would also write it up in a blog post.

1 Overview

I am going to describe two design patterns, “object algebras” and “finally tagless interpreters”, developed by the object-oriented and functional communities to solve the expression problem. I hope my presentation will serve as a hands-on, down-to-earth introduction to both patterns, but more importantly I am going to argue that they are really implementations of the same idea.

The plan is as follows. First, we are going to recall what the expression problem is about. Then, I will introduce object algebras and we will see how they help to solve the expression problem. Finally, we will see that the Java code for object algebras admits a pretty natural translation to Haskell, which will lead us to the finally tagless encoding.

2 Expression Problem

We begin by discussing the expression problem. To this end, we will use a classic example known as “Hutton’s Razor”. It is a language of simple arithmetic expressions built up from integer literals using the addition operation +. In Haskell, expressions in this language are described by the following algebraic data type Exp:

data Exp = Lit Int | Add Exp Exp

For example, the expression (1 + (2 + 3)) is represented in Haskell as follows:

e1 = Add (Lit 1)
         (Add (Lit 2)
              (Lit 3))

The function eval takes an expression as an argument and returns its integer value:

eval :: Exp -> Int
eval (Lit n)   = n
eval (Add x y) = eval x + eval y

This is a minimal setup that we will try to extend in various ways.

Suppose that we want to add a function view that converts an expression to a string, its textual representation. This is fairly straightforward:

view :: Exp -> String
view (Lit n)   = show n
view (Add x y) = "(" ++ view x ++ " + " ++ view y ++ ")"

If, however, we decide to extend our language with the multiplication operation *, we have to introduce a new constructor Mul in the definition of the data type Exp, which entails adding a new case to every function defined on Exp.

Java, as well as many other object-oriented languages, doesn’t have algebraic data types and the standard object-oriented approach to this problem is different.

To wit, one defines an interface (or an abstract class) Exp of expressions, which is implemented (or subclassed) by concrete classes Lit and Add representing different types of expressions:

interface Exp {
    int eval();
}

class Lit implements Exp {
    int n;
    int eval() { return n; }
}

class Add implements Exp {
    Exp x, y;
    int eval() { return x.eval() + y.eval(); }
}

For the sake of brevity, from now on I will omit access modifiers and constructor definitions in Java code.

The operation eval is defined as a method of the interface Exp, and the concrete classes Lit and Add implement this method.

In order to add the multiplication operation, it suffices to define a new class Mul that implements the interface Exp:

class Mul implements Exp {
    Exp x, y;
    int eval() { return x.eval() * y.eval(); }
}

Note that we don’t need to touch the classes Lit and Add.

However, it is rather awkward to add new operations as we need to add new methods to the interface Exp, which is going to break all its concrete implementations.

We see that the functional and object-oriented paradigms are in some sense dual to each other. There is a standard way to visualize this duality.

Consider a rectangular table whose columns are labeled by operations and whose rows are labeled by expression types, and the cell at the intersection of given column and row contains a definition of the corresponding operation for the corresponding expression type.

2014-06-19-23:29:34.png

In the functional approach, the definitions in the table are grouped by columns. A column corresponding to an operation represents a function containing definitions of that operation for all types of expressions. The table is naturally extensible in the dimension of operations: adding a new operation is easy and corresponds to adding a new function, i.e., a new column. Extensibility in the dimension of expression types is problematic: adding a new expression type is hard as it corresponds to adding a new row crossing all columns, and every intersection is a new case that has to be added to the corresponding function definition.

2014-06-19-23:29:42.png

In the object-oriented approach, the definitions are grouped by rows. A row corresponding to a given expression type represents a class containing definitions of all operations for that type. One of the promises of the object-oriented programming is painless extensibility in the direction of expression types. Indeed, adding a new expression type is easy in the object-oriented approach and corresponds to adding a new class, i.e., a new row to the table. Adding a new operation is hard as it corresponds to adding a new column crossing all rows, and every intersection is a new method definition that has to be added to each of the existing classes.

2014-06-19-23:29:49.png

Can we achieve extensibility in both dimensions without loss of static guarantees and without having to modify existing code? That is the essence of the expression problem.

3 Object Algebras

Object algebras is a relatively new object-oriented design pattern for solving the expression problem. It has been introduced in the paper “Extensiblity for the Masses”. What makes it attractive is that it is a relatively lightweight pattern that doesn’t require fancy language extensions. In particular, it can be easily implemented in mainstream languages like Java and C#. Let us see how.

We begin again by considering the language containing only literals and addition. It is going to be described by an interface

interface ExpAlg<T> {
    T lit(int n);
    T add(T x, T y);
}

If you are familiar with GoF design patterns, you may recognize this as the Visitor interface. However, in this approach this interface is going to play a completely different role. Namely, in terms of GoF design patterns ExpAlg<T> is an interface of an abstract factory for creating expressions.

For example, the expression (1 + (2 + 3)) is represented as follows:

<T> T e1(ExpAlg<T> f) {
    return f.add(
        f.lit(1),
        f.add(
            f.lit(2),
            f.add(3)));
}

That is, it is represented as a function taking as an argument an object f implementing the interface ExpAlg<T> (i.e., a concrete factory) and returning a value of type T. The body of the function simply calls the methods lit and add of the factory f in a suitable order with suitable arguments.

This representation allows us to vary the concrete factory f thus interpreting the expression e1 in different ways.

Let us see how this works in the case of evaluation of expressions.

First of all, we introduce an interface Eval of “objects that can be evaluated”:

interface Eval { int eval(); }

Next we define a concrete factory EvalExp, which is going to manufacture expression that can be evaluated:

class EvalExp implements ExpAlg<Eval> {
    Eval lit(final int n) {
        return new Eval() {
            int eval() {
                return n;
            }
        };
    }
    Eval add(final Eval x, final Eval y) {
         return new Eval() {
             int eval() {
                 return x.eval() + y.eval();
             }
         };
    }
}

We can now pass an instance of the factory EvalExp into the expression e1 and get back an object that has a method eval. We can call this method and convince ourselves that it returns the right value:

int v1 = e1(new EvalExp()).eval();

We shall soon see another example of defining operations, but first let us think how we could add multiplication to our language.

We could add a new method mul to the interface ExpAlg<T>, but then the implementation of the concrete factory EvalExp would require changes, which is precisely what we would like to avoid.

Instead, we introduce a new interface MulAlg<T> that extends the interface ExpAlg<T> and adds a new method mul to it:

interface MulAlg<T> extends ExpAlg<T> {
    T mul(T x, T y);
}

Expressions containing multiplication are now going to be represented as functions taking as an argument objects implementing the extended interface MulAlg<T>. For example, the expression (4 * (5 + 6)) will be represented as follows:

<T> T e2(MulAlg<T> f) {
    return f.mul(
        f.lit(4),
        f.add(
            f.lit(5),
            f.lit(6)));
}

To extend the implementation of evaluation of expressions to expressions containing multiplication we define a new concrete factory EvalMul that implements the interface MulAlg<Eval> and inherits from the factory EvalExp implementations of the methods lit and add:

class EvalMul extends EvalExp implements MulAlg<Eval> {
    Eval mul(final Eval x, final Eval y) {
        return new Eval() {
            int eval() {
                return x.eval() * y.eval();
            }
        };
    }
}

We can now pass an instance of the factory EvalMul into the expression e2, get back an object that can be evaluated, and compute its value by calling the eval method:

int v2 = e2(new EvalMul()).eval();

Note that we are not touching any existing code: we are defining new interfaces and classes and use inheritance to avoid duplication.

Let us consider an example of adding a new operation and convince ourselves that it doesn’t require touching existing code either.

We introduce an interface View of “object that can be converted to a string”:

interface View { String view(); }

We then define concrete factories ViewExp and ViewMul implementing the interfaces ExpAlg<View> and MulAlg<View>:

class ViewExp implements ExpAlg<View> {
    View lit(final int n) {
        return new View() {
            String view() {
                return Integer.toString(n);
            }
        };
    }
    View add(final View x, final View y) {
        return new View() {
            String view() {
                return "(" + x.view() + " + " + y.view() + ")";
            }
        };
    }
}

class ViewMul extends ViewExp implements MulAlg<View> {
    View mul(final View x, final View y) {
        return new View() {
            String view() {
                return "(" + x.view() + " * " + y.view() + ")";
            }
        };
    }
}

Notice that we have a choice: to define conversion to a string only for expression containing literals and addition or also for expressions containing multiplication. In the latter case, the class ViewMul inherits from ViewExp the methods responsible for conversion of literals and addition and adds a definition of conversion of multiplication, thereby implementing the interface MulAlg<View>.

We can now pass these factories into the expressions e1 and e2, get back objects supporting the method view, and convince ourselves that this method returns the right values:

String s1 = e1(new ViewExp()).view();
String s2 = e2(new ViewMul()).view();

A few words about the terminology. The abstract factory interface ExpAlg<T> is called an object algebra interface and concrete factories implementing it are called object algebras. The terminology and the approach in general are inspired by abstract algebra, in which algebraic structures (or simply algebras) are described by signatures. A signature acts as an interface that specifies the types of operations defined on the underlying set of an algebraic structure, and an algebra provides concrete definitions of those operations and is similar to a class implementing an interface.

4 Finally Tagless

We now translate the above example from Java to Haskell and arrive at essentially the finally tagless approach.

This approach has been introduced in the paper “Finally Tagless, Partially Evaluated”. It enjoys many interesting properties, but we are only going to focus on how it helps us to solve the expression problem.

First, we consider how expressions are represented in the finally tagless approach.

The Haskell counterpart of the interface ExpAlg<T> in Java is a typeclass ExpAlg:

class ExpAlg t where
    lit :: Int -> t
    add :: t -> t -> t

The expression (1 + (2 + 3)) is represented as follows:

e1 = add (lit 1)
         (add (lit 2)
              (lit 3))

and its inferred type ExpAlg t => t can be read as a statement that the term e1 can be given any type t as long as t is an instance of the typeclass ExpAlg, i.e., supports the operations lit and add.

What is remarkable is that at the stage of desugaring the typeclass ExpAlg is translated into a record type with fields lit and add, instance declarations of ExpAlg are translated into record values of this type, and the expression e1 is translated into a function taking as an argument a record value corresponding to an instance declaration of ExpAlg for a given type t. This record value is a direct analog of the concrete factory f in Java.

The difference between Haskell and Java is that in Haskell for every pair of a typeclass and a type there can be at most one instance declaration, and once this instance declaration is defined, it remains implicitly in the scope and is passed around automatically by the compiler. In Java, there can be many implementations of a given generic interface instantiated with a given type, and these implementation are entities in the program that have names and have to be passed around explicitly.

To understand how operations are defined in the finally tagless approach, let us look again at evaluation of expressions. We are going to simplify the Java implementation first to make translation to Haskell more obvious.

Namely, we observe that we don’t have to make Eval an interface, although conceptually this is cleaner. After all, if all we know about an object is that it implements the interface Eval, then all we can do with this object is to call the method eval. Therefore, for an external observer such an object is practically indistinguishable from an object of a class Eval with a public field eval of type int. With this change, the definition of the factory EvalExp becomes:

class Eval { public int eval; }

class EvalExp implements ExpAlg<Eval> {
    Eval lit(int n) {
        return Eval(n);
    }
    Eval add(Eval x, Eval y) {
        return Eval(x.eval + y.eval);
    }
}

To evaluate the expression e1 we pass it an instance of the class EvalExp, which produces an object of the class Eval, which we can ask the value of its field eval:

int v1 = e1(new EvalExp()).eval;

The class Eval is a wrapper around the type int. In Haskell, we introduce a wrapper Eval around the type Int implemented, say, as a newtype:

newtype Eval = Eval { eval :: Int }

The definition of the class EvalExp translates into an instance declaration of the typeclass ExpAlg for the type Eval:

instance ExpAlg Eval where
    lit n   = Eval n
    add x y = Eval $ eval x + eval y

To evaluate the expression e1 we restrict its type ExpAlg t => t to Eval, which is similar to applying the expression e1 to the concrete factory EvalExp in Java, and then the obtained value of type Eval can be queried the value of the field eval:

v1 = eval (e1 :: Eval)

In fact, we don’t have to explicitly restrict the type of e1 to Eval as the compiler will infer this automatically from the fact that e1 is passed as an argument to the function eval :: Eval -> Int.

In Java, extension of the language with the multiplication operation corresponded to extension of the interface ExpAlg<T> to MulAlg<T> with the method mul. In Haskell, this translates into a definition of a new typeclass MulAlg that is a subclass of ExpAlg:

class ExpAlg t => MulAlg t where
    mul :: t -> t -> t

Expressions containing multiplication will now have the type MulAlg t => t, for example:

e2 = mul (lit 4)
         (add (lit 5)
              (lit 6))

To extend the definition of the operation eval, we perform the same transformation as above: we replace the interface Eval with a class Eval with a public field of type int. Then the definition of the class EvalMul implementing the interface MulAlg<Eval> translates in Haskell into a definition of an instance declaration of the typeclass MulAlg for the type Eval:

instance MulAlg Eval where
    mul x y = Eval $ eval x * eval y

Finally, we consider an example of adding new operations in the finally tagless approach. Let us implement conversion to string.

We transform the object algebras implementation of conversion in Java similarly to what we did to the implementation of evaluation: we replace the interface View with a class View with a public field view of the type String. Then in Haskell this all is going to turn into an instance declaration of the typeclass ExpAlg for a type View that is a newtype wrapper around String:

newtype View = View { view :: String }

instance ExpAlg View where
    lit n   = View $ show n
    add x y = View $ "(" ++ view x ++ " + " ++ view y ++ ")"

Extension of the operation view to expressions containing multiplication is done similarly and amounts to defining an instance declaration of the typeclass MulAlg for the type View:

instance MulAlg View where
    mul x y = View $ "(" ++ view x ++ " * " ++ view y ++ ")"

The encoding we have arrived at is almost identical to that proposed in “Finally Tagless, Partially Evaluated”. The difference is that in Haskell there is no need to make MulAlg a subclass of ExpAlg: we can define it as a separate typeclass:

class MulAlg t where
    mul :: t -> t -> t

and then expressions containing multiplication will be represented by values of the type (ExpAlg t, MulAlg t) => t. This gives us more flexibility because we can represent different language extensions as separate typeclasses that can be arbitrarily combined.

CPSing anamorphisms

At the last meeting of the London Haskell User Group, my boss Peter Marks gave a talk in which he presented definitions of catamorphisms and anamorphisms such that the composition of an anamorphism followed by a catamorphism (also known as a hylomorphism) doesn’t build up an intermediate data structure.

Even though the definitions are pretty neat, the idea is not new and boils down to CPSing the definition of anamorphism. In this blog post, I’d like to illustrate on the example of lists how the CPSed definition can be derived from the standard one with the help of the Yoneda embedding (and thereby also illustrate the idea that CPS and the Yoneda embedding are the same thing).

We begin by recalling the definitions of catamorphism and anamorphism in the case of lists. We represent lists of elements of type a as the fixed point of the functor ListF a:

data ListF a b = Nil | Cons a b

instance Functor (ListF a) where
  fmap _ Nil        = Nil
  fmap f (Cons a b) = Cons a (f b)

newtype Fix f = Fix { unFix :: f (Fix f) }

type List a = Fix (ListF a)

Catamorphisms and anamorphisms are then defined as follows:

cata :: (ListF a b -> b) -> List a -> b
cata alg = alg . fmap (cata alg) . unFix

ana :: (b -> ListF a b) -> b -> List a
ana coalg = Fix . fmap (ana coalg) . coalg

In his talk, Peter used a different representation of lists, which he called a list abstraction:

data ListAbs e v a = ListAbs (e -> v -> a) a

instance Functor (ListAbs e v) where
  fmap f (ListAbs c n) = ListAbs (\e v -> f (c e v)) (f n)

nil :: ListAbs e v a -> a
nil (ListAbs _ n) = n

cons :: e -> v -> ListAbs e v a -> a
cons e v (ListAbs c _) = c e v

and defined catamorphisms and anamorphisms as follows:

type ListAlg e a = ListAbs e a a

cataList :: ListAlg e a -> [e] -> a
cataList alg = go
  where
    go (x : xs) = cons x (go xs) alg
    go []       = nil            alg

type ListCoalg v e = forall a. v -> ListAbs e v a -> a

anaList :: ListCoalg v e -> v -> [e]
anaList f v = f v x
  where
    x = ListAbs (\e v -> e : f v x) []

The type ListF a b -> c is clearly isomorphic to the type ListAbs a b c:

p :: (ListF a b -> c) -> ListAbs a b c
p alg = ListAbs (\e v -> alg (Cons e v)) (alg Nil)

q :: ListAbs a b c -> (ListF a b -> c)
q abs (Cons e v) = cons e v abs
q abs Nil        = nil      abs

Therefore, the type ListF a b -> b is isomorphic to ListAbs a b b, which by definition is equal to ListAlg a b. This isomorphism allows us to translate the function cata into a function

cata' :: ListAlg a b -> List a -> b
cata' = cata . q

Substituting the definitions of cata and q and simplifying, we obtain:

cata' alg = go
  where
    go (Cons x xs) = cons x (go xs) alg
    go Nil         = nil            alg

In other words, we arrive at the definition of cataList, up to the isomorphism between List a and a.

To relate ana and anaList in a similar way, we shall use the Yoneda Lemma, which claims that the type forall b. (a -> b) -> f b of natural transformations from a representable functor (->) a to an arbitrary functor f is isomorphic to f a. The isomorphism is given by a pair of mutually inverse functions:

y :: (forall b. (a -> b) -> f b) -> f a
y f = f id

z :: Functor f => f a -> (forall b. (a -> b) -> f b)
z = flip fmap

In particular, the type forall c. (b -> c) -> (a -> c) of natural transformations from a representable functor (->) b to a representable functor (->) a is isomorphic to the type a -> b of functions from a to b.

By the Yoneda Lemma, the type b -> ListF a b is isomorphic to the type forall c. (ListF a b -> c) -> b -> c. The type ListF a b -> c is isomorphic to ListAbs a b c. Therefore, b -> ListF a b is isomorphic to forall c. ListAbs a b c -> b -> c, which up to the order of arguments is ListCoalg a b. This isomorphism allows us to translate the function ana into a function

ana' :: (forall c. ListAbs a b c -> b -> c) -> a -> List b
ana' = ana . y . (. p)

Expanding the definitions, we obtain:

ana' f = ana ((f . p) id)
       = ana (f (p id))
       = ana (f (ListAbs Cons Nil))

This doesn’t resemble anaList yet. Here is a trick. Let us denote the function f (ListAbs Cons Nil) :: b -> ListF a b by g:

ana' f = ana g
  where
    g = f (ListAbs Cons Nil)

and do a bit of equational reasoning:

ana g

  -- definition of ana
  = Fix . fmap (ana g) . g

  -- definition of g
  = Fix . fmap (ana g) . f (ListAbs Cons Nil)

  -- definition of fmap for (->) b
  = fmap (Fix . fmap (ana g)) . f $ ListAbs Cons Nil

  -- naturality of f
  = f . fmap (Fix . fmap (ana g)) $ ListAbs Cons Nil

  -- definition of fmap for ListAbs and ListF
  = f $ ListAbs (\e v -> Fix (Cons e (ana g v))) (Fix Nil)

The penultimate equality is by the naturality of f, which is a consequence of the free theorem associated with the polymorphic type forall c. ListAbs a b c -> b -> c. We conclude that

ana' f = h
  where
    h = f $ ListAbs (\e v -> Fix (Cons e (h v))) (Fix Nil)

and you can easily convince yourself that this definition is equivalent to that of anaList.

Ariadne plugin for Emacs

I am pleased to announce ariadne-el, an Emacs interface to Ariadne.

Ariadne is a new tool written by Roman Cheplyaka that provides the go-to-definition functionality for Haskell code. It is designed as a server that responds to queries from IDEs and text editor plugins transmitted via the BERT-RPC protocol over TCP. The server is implemented in Haskell using Roman’s haskell-names name resolution library. It is fully Haskell-aware and can properly locate prefixed names (like T.head) and locally bound names, which makes it quite a bit smarter than TAGS. Ariadne is in an early development stage, and as of v0.1.2 it can only find definitions in the current file, but I am sure this and other limitations will be lifted in future versions.

ariadne-el allows you to communicate with Ariadne from Emacs. It uses my own bert-el, BERT serialization library for Emacs, for encoding/decoding BERT-RPC protocol messages. I have to admit that I never worked with networking facilities of Emacs before, so I shamelessly borrowed the design of low-level networking from SLIME.

If you are a Haskell programmer and an Emacs user, I encourage you to give Ariadne and ariadne-el a try. The installation is pretty straightforward:

  1. Install Ariadne and run the ariadne-server executable.
  2. Install bert-el and ariadne-el by cloning the Git repositories and adding their paths to load-path. If you use Marmalade, you can also obtain both packages by running M-x package-install ariadne.
  3. Add the following lines to your .emacs file:
(require 'ariadne)
(add-hook 'haskell-mode-hook
          (lambda ()
            (define-key haskell-mode-map "\C-cd" 'ariadne-goto-definition)))

Now try to open any Haskell source file, put the cursor on any name, and press C-c d.

Co-Yoneda Lemma

The purpose of this humble blog post is to understand this insightful comment by Edward Kmett to the question “Is operational really isomorphic to a free monad?” posted at StackOverflow.

Consider the following data type definition in Haskell, which can also be found in Edward Kmett’s package kan-extensions:

data CoYoneda f a = forall b. CoYoneda (b -> a) (f b)

In Haskell, a forall before the data constructor means existential rather than universal quantification, so that CoYoneda f a is isomorphic to the type \exists b. (b\to a, f\, b). The forall hints at the polymorphic type of the constructor

CoYoneda :: forall b. (b -> a) -> f b -> CoYoneda f a

Every polymorphic function satisfies a theorem that can be extracted by reading its type as a relation. This is referred to as parametricity. The reader is referred to Wadler’s classic paper “Theorems for free!” for more details.

Suppose that f is a functor. Then parametricity implies the following free theorem for the function CoYoneda: if g :: b -> a and y :: f c, then CoYoneda (g . h) y == CoYoneda g (fmap h y) for any h :: c -> b.

Consider the following two functions:

phi :: f a -> CoYoneda f a
phi x = CoYoneda id x

psi :: Functor f => CoYoneda f a -> f a
psi (CoYoneda g y) = fmap g y

Then obviously psi . phi == id:

(psi . phi) x == fmap id x == x

The above free theorem implies also that phi . psi == id:

(phi . psi) (CoYoneda g y) == CoYoneda id (fmap g y) == CoYoneda g y

We conclude that CoYoneda f is isomorphic f. This is a Haskell version of the co-Yoneda Lemma.

Appendix

The only thing that is missing in Wadler’s paper is how to associate a relation with the type f b. Wadler only explains this in the case of the list functor. We assume that the Haskell functor f is modeled by a functor F on the category of sets. Suppose that \mathcal{B}: B \Leftrightarrow B' is a relation associated with the type b. That is, \mathcal{B} is a subset of B\times B'. Let \iota: \mathcal{B}\hookrightarrow B\times B' denote the inclusion map. We define the relation associated with the type f b as the image of the composite \langle F\pi, F\pi'\rangle\circ F\iota: F(\mathcal{B})\to FB\times FB', where \pi: B\times B'\to B and \pi': B\times B'\to B' are the canonical projections.

For example, if f = [], so that F=(-)^* is the list functor, then the above definition says that two lists xs \in B^* and ys \in B'^* are related if there exists a list zs\in \mathcal{B}^* \subset (B\times B')^* such that

\displaystyle  \begin{array}{rcl} map\; fst\; zs & \equiv & xs,\\ map\; snd\; zs & \equiv & ys. \end{array}

In other words, xs and ys are related if they have the same length and corresponding elements are related, which coincides with Wadler’s definition.

It follows from this definition that if \mathcal{B} is the relation associated with a function g: B\to B', then F(\mathcal{B}) is the relation associated with the function Fg. Indeed, if \mathcal{B} is the graph of g, i.e., the image of the map \langle\mathrm{id}_B, g\rangle: B\to B\times B', then F(\mathcal{B}) is the image of the composite \langle F\pi, F\pi'\rangle\circ F(\langle\mathrm{id}_B, g\rangle), which coincides with \langle F\mathrm{id}_B, Fg\rangle = \langle \mathrm{id}_{FB}, Fg\rangle. Therefore, F(\mathcal{B}) coincides with the graph of Fg.

Category Theory Helpline

One of my colleagues, Ben Moseley, has recently asked an interesting category theoretic question. I’ve decided to write up my answer in a blog post.

Preliminaries

This blog post is a literate Haskell program. First, we are going to enable a few language extensions and import a few modules:

> {-# LANGUAGE FlexibleInstances     #-}
> {-# LANGUAGE MultiParamTypeClasses #-}
> {-# LANGUAGE TupleSections         #-}
> {-# LANGUAGE TypeOperators         #-}
> import Control.Arrow
> import Data.Traversable
> import Prelude hiding (sequence)

Next, we are going to recall the definitions of the fixed point of a functor f, the free monad generated by f, and the cofree comonad generated by f:

> newtype Fix f = Fix { unFix :: f (Fix f) }
> data Free f a = Pure a | Free (f (Free f a))
> data Cofree f a = a :< f (Cofree f a)

All three definitions are fairly standard. The last two can be found, for example, in Edward Kmett’s package free.

Fix f is both an initial algebra and final coalgebra for f: the f-algebra structure is given by Fix :: f (Fix f) -> Fix f, and the f-coalgebra structure is given by unFix :: Fix f -> f (Fix f).

The initiality means that for each f-algebra alg :: f a -> a there exists a unique f-algebra morphism cata alg :: Fix f -> a, called catamorphism. cata alg being an f-algebra morphism means that the following equation holds:

cata alg . Fix == alg . fmap (cata alg)

from which we can read off the definition of cata alg:

> cata :: Functor f => (f a -> a) -> Fix f -> a
> cata alg = alg . fmap (cata alg) . unFix

Similarly, the finality means that for each f-coalgebra coalg :: a -> f a there exists a unique f-coalgebra morphism ana coalg :: a -> Fix f, called anamorphism. ana coalg being an f-coalgebra morphism means that the following equation holds:

unFix . ana coalg == fmap (ana coalg) . coalg

from which we can read off the definition of ana coalg:

> ana :: Functor f => (a -> f a) -> a -> Fix f
> ana coalg = Fix . fmap (ana coalg) . coalg

The original version of this blog post fell into the fallacy of claiming that every function from resp. to Fix f could be expressed a catamorphism resp. an anamorphism. This is obviously not the case. “When is a function a fold or an unfold?” by Jeremy Gibbons, Graham Hutton, and Thorsten Altenkirch gives necessary and sufficient conditions for when a set-theoretic function can be written as a catamorphism resp. an anamorphism (thanks to Ben for pointing this out). The fallacy doesn’t affect the calculations below. However, they become somewhat less natural than I originally thought.

The functions cata and ana are only two examples of the plethora of recursion schemes associated with the type Fix f. The types Free f a and Cofree f a also have associated recursion schemes:

> fold :: Functor f => (Either a (f b) -> b) -> Free f a -> b
> fold f m = f $ case m of
>   Pure x -> Left x
>   Free n -> Right $ fmap (fold f) n

> unfold :: Functor f => (b -> (a, f b)) -> b -> Cofree f a
> unfold f c = case f c of
>   (x, d) -> x :< fmap (unfold f) d

The question

Ben wrote the following functions:

> cofree2fix :: Functor f => Cofree f a -> Fix f
> cofree2fix = ana (\(_ :< fx) -> fx)

> fix2free :: Functor f => Fix f -> Free a
> fix2free = cata Free

> fix2cofree :: Functor f => (a, Fix f) -> Cofree f a
> fix2cofree (a, fx) = unfold (const a &&& unFix) fx

> free2fix :: Traversable f => Free f a -> Either a (Fix f)
> free2fix = fold (Left ||| (right Fix . sequence))

Here

(&&&) :: (a -> b) -> (a -> c) -> a -> (b, c)

and

(|||) :: (b -> d) -> (c -> d) -> Either b c -> d

are the standard mediating morphisms for products and sums. They are imported from Control.Arrow. The function sequence is the generalized version from Data.Traversable.

Ben then wrote:

I’m puzzled by the last one… It all typechecks, but there’s a clear asymmetry between fix2cofree and free2fix which worries me…

The rest of this blog post is an attempt to explain what is going on here.

It’s fixed points all the way down

We first observe that Free f a is isomorphic to the fixed point of the functor Either a :. f, and Cofree f a is isomorphic to the fixed point of the functor (,) a :. f, where :. denotes functor composition (thanks to Roman Cheplyaka for suggesting this notation):

> newtype (f :. g) a = C { unC :: f (g a) }

> instance (Functor f, Functor g) => Functor (f :. g) where
>   fmap t (C x) = C $ fmap (fmap t) x

We can encode isomorphisms between types by means of a typeclass, Iso:

> class Iso a b where
>   iso :: a -> b
>   inv :: b -> a

The isomorphisms between Free f a and Fix (Either a :. f), and between Cofree f a and Fix ((,) a :. f) are then defined as follows:

> instance Functor f => Iso (Free f a) (Fix (Either a :. f)) where
>   iso = fold (Fix . C)
>   inv = cata (either Pure Free . unC)

> instance Functor f => Iso (Cofree f a) (Fix ((,) a :. f)) where
>   iso = ana (\(x :< fx) -> C (x, fx))
>   inv = unfold (unC . unFix)

Under this identification, the recursion schemes fold and unfold identify with cata and ana respectively.

We shall show now that the functions cofree2fix, fix2free, fix2cofree, and free2fix can naturally be expressed as cata- resp. anamorphisms of some naturally arising algebras resp. coalgebras over the functor f.

cofree2fix

The finality of Fix f allows us to manufacture a function

Fix ((,) a :. f) -> Fix f

out of a function

Fix ((,) a :. f) -> f (Fix ((,) a :. f))

There is an obvious natural representative of the latter type, namely the composite snd . unC . unFix, leading to the following definition:

> cofree2fix' :: Functor f => Fix ((,) a :. f) -> Fix f
> cofree2fix' = ana (snd . unC . unFix)

fix2free

The initiality of Fix f implies that each function

f (Fix (Either a :. f)) -> Fix (Either a :. f)

gives rise to a function

Fix f -> Fix (Either a :. f)

There is an obvious natural representative of the former type, namely the composite Fix . C . Right, which corresponds to the constructor Free under the isomorphism between Free f a and Fix (Either a :. f). Therefore, we can define

> fix2free' :: Functor f => Fix f -> Fix (Either a :. f)
> fix2free' = cata (Fix . C . Right)

fix2cofree

By the finality of Fix ((,) a :. f), we obtain a function

(a, Fix f) -> Fix ((,) a :. f)

from each function

(a, Fix f) -> ((,) a :. f) (a, Fix f)

The type ((,) a :. f) (a, Fix f) is isomorphic to (a, f (a, Fix f)) by means of the constructor C. Therefore, we need to define a function of type

(a, Fix f) -> (a, f (a, Fix f))

By the universal property of product, any such function is necessarily of the form u &&& v for uniquely determined

u :: (a, Fix f) -> a

and

v :: (a, Fix f) -> f (a, Fix f)

There is an obvious candidate for u: the function fst (and, in fact, because we are after a polymorphic u, this is the only choice we have). Is there also a natural choice of v? The answer is yes, and the fundamental reason for that is that Haskell functors are strong. That is, an arbitrary functor f admits the following natural transformation, which in category theory is called the right tensorial strength:

> strength :: Functor f => (a, f b) -> f (a, b)
> strength (x, fy) = fmap (x,) fy

This allows us to define

> fix2cofree' :: Functor f => (a, Fix f) -> Fix ((,) a :. f)
> fix2cofree' = ana (C . (fst &&& strength . second unFix))

free2fix

The initiality of Fix (Either a :. f) implies that each function

(Either a :. f) (Either a (Fix f)) -> Either a (Fix f)

gives rise to a function

Fix (Either a :. f) -> Either a (Fix f)

Therefore, in order to construct a function of the latter type it suffices to come up with a function of the former type.

The type (Either a :. f) (Either a (Fix f)) is isomorphic to Either a (f (Either a (Fix f))) by means of C, therefore we are after a function of type

Either a (f (Either a (Fix f))) -> Either a (Fix f)

By the universal property of sum, any such function is necessarily of the form p ||| q with uniquely determined

p :: a -> Either a (Fix f)

and

q :: f (Either a (Fix f)) -> Either a (Fix f)

We have an obvious candidate for p: the function Left (and because p has to be polymorphic, this is the only choice we have). By analogy with the previous case, we might expect that there always exists a natural transformation of type

Functor f => f (Either a b) -> Either a (f b)

Alas, this is not the case. There is, however, a large class of functors that do admit a natural transformation of this type: Traversable functors, and the required function is sequence.

Therefore, for Traversable functors we can define

> free2fix' :: Traversable f => Fix (Either a :. f) -> Either a (Fix f)
> free2fix' = cata ((Left ||| (right Fix . sequence)) . unC)

Summary

The short answer to Ben’s question is this: fix2cofree works for any functor because all functors in Haskell are strong: for each functor f, there is a natural transformation strength :: Functor f => (a, f b) -> f (a, b), called the right tensorial strength, subject to coherence axioms.

The fundamental fact underlying the existence of strength is that, categorically, every Haskell functor is enriched: its action on morphisms is given not merely as a map from the set of functions a -> b to the set of functions f a -> f b, but as a morphism (in the category Hask) from the internal object of morphisms from a to b, the type a -> b, to the internal object of morphisms from f a to f b, the type f a -> f b. This morphism is fmap, of course. There is a one-to-one correspondence between enrichments of a functor F from a cartesian closed category \mathcal{C} to itself and right tensorial strengths on F. For more information on strong functors, see this nLab web page, or refer to the original paper by Anders Kock.

On the other hand, not every Haskell functor f admits a dual natural transformation Functor f => f (Either a b) -> Either a (f b). Traversable functors do, but the former condition is weaker, I think.

Update: User sclv on reddit suggested that (up to the order of summands) the function f (Either a b) -> Either a (f b) is called costrength in Edward Kmett’s post.

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]

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.

My Setup

On and off, I read interviews about people’s setups at usesthis.com. I particularly enjoyed the interviews of Andrew Huang, Benjamin Mako Hill, Phil Hagelberg, and Russ Cox. For amusement, I’ve decided to reflect on my own computer usage and compile a description of my setup.

Who are you, and what do you do?

I’ve been trained as a mathematician, but over the last a few years I’ve become more interested in programming, particularly functional and mathematically structured programming, and in theoretical computer science, most prominently programming languages theory and compilers. I’ve spent the last two years as a PhD student at NUIM thinking primarily about design and implementation of a novel functional programming language for scientific computing with support for automatic differentiation. I’m beginning a new job in industry in September. I’m going to work on a large commercial product written in Haskell, and I’m looking forward to it.

What hardware do you use?

My only machine is a Lenovo ThinkPad SL 500 laptop running a 2.00 GHz Intel Core 2 Duo T5870, 3Gb of RAM. It is a two years old machine, and its age is starting to show. Still, at the moment of purchase, it was probably the best laptop with no Windows pre-installed one could buy in Ukraine.

I have recently acquired a Kindle 4, non-touch edition, which I use for reading fiction. I have many papers and books in PDF format that I’d also like to read away from my computer, but I haven’t found an affordable device suitable for that.

And what software?

I run Ubuntu 10.04 (64 bit) on my laptop. I’ve been reluctant to upgrade because upgrades are know to break things and because this particular version has been working pretty well for me. When I want the latest version of some software, I build it from source. The following is a list of applications I use on a daily basis, in the loose order of descending importance.

Emacs

When I’m working (and often also when I’m not), I’m spending most of my time in Emacs. Compared to other Emacs users, the severity of my addiction to Emacs is modest: Emacs is my text editor, organizer, file manager, and shell. I also have an Emacs interface to Google Translate, which I use as a poor man’s dictionary.

Emacs is one of few pieces of software that I build from source.

Browser

The second most used application is the browser. I’ve tried a few different ones: Firefox, Google Chrome, Opera, Chromium, Conkeror. I’m not happy with any of them. The browser I’ve enjoyed most was Conkeror, which is to browsers what Emacs is to editors. However, I had to abandon it because occasionally it would blow up and consume all available RAM and half of available swap, and because it reliably crashes when the laptop is awaken from suspend mode.

Right now I’m using Chromium, for no particular reason. I try to keep the number of open tabs low, otherwise Chromium as well as Google Chrome become a memory drain on my machine.

xmonad

My window manager is xmonad. A lot has been said in praise of tiling window managers, in particular, about how they make you more productive than conventional desktop environments. I’m not sure I buy that.

I started to use xmonad mainly out of curiosity and because of the hype surrounding it, and also because it is implemented in Haskell, so getting into it was an opportunity to learn some Haskell. Later, when I was interning at Streamtech, a webdev startup from The Hague, I enjoyed the ability to easily tile the browser and a bunch of terminal windows on my screen. Still later, when I was doing PhD at NUIM, I discovered what a joy it was to use xmonad with a dual head setup.

However, right now I’m using it only on my laptop, and I’m not sure if it’s a such a big win there. Unlike other xmonad users, I’m not a terminal junkie. In fact, I don’t run terminals outside Emacs at all, except for the case when I want to ssh to a remote host. I also rarely tile my windows because I don’t have enough screen real estate for that. I run pretty much every application fullscreen all the time, with only a few exceptions. I do like the idea of having dedicated workspaces for different activities, but I rarely use more than 3 simultaneously.

Running pure xmonad on a laptop is also somewhat awkward because one has to think about many things DEs take care for you, for example, having different keyboard layouts and a layout indicator, having a system tray with the network manager applet and battery status, mounting devices etc. That’s why I run xmonad as a drop-in replacement for Metacity in Gnome. This way I get the best of both worlds.

However, I don’t like the direction in which Gnome is heading, and I’ve never liked KDE, so most likely I’m stuck with xmonad, and should I need to reinstall my system, I’ll probably go for some combination of xmonad and Gnome tools.

git

I use git as my version control system. I almost never use the command line interface. Instead, for frequent tasks like stage, commit, push, pull, create/change/delete a branch, and view diff/log I use magit, an awesome Emacs interface to git. For more complex git surgery there is gitk. I also use git-gui when I want to split a large commit into smaller chunks.

LaTeX

I use LaTeX for typesetting my papers. AUCTeX mode for Emacs is superb. I’m surprised to see people using both LaTeX and Emacs and not using AUCTeX, but rather invoking LaTeX, BibTeX, dvitops and other tools from a shell, either manually or using a makefile. Those people don’t know what they are missing.

Evince

I use Evince to view my PDF, PostScript, and DJVU files. Evince is not particularly configurable, and I hate that I have to change settings for every new open file, but it does its job sufficiently well that I don’t have an urge to switch to something better.

My only gripe about Evince is that its DVI viewer doesn’t support forward/inverse search.

xdvi

That’s why I use xdvi when working with LaTeX. xdvi is old and ugly as hell, but it is fast and its forward/inverse search are killer features for me. Making inverse search work with Emacs used to be a PITA, but AUCTeX takes care of that, too.

Skype

I use Skype for VOIP. I like how much cleaner Skype looks on Linux compared to the bloated interface on Windows. For example, I was unable to share screen using Skype on my wife’s Windows laptop.

Calibre

As I mentioned above, I bought myself a Kindle 4 recently. Because I didn’t want to deal with Amazon, I needed an application to conveniently upload e-books to my device and to occasionally convert between different formats. Calibre performs these tasks seamlessly.

Compilers and Interpreters

I like to explore programming languages, and I have a few compilers and interpreters installed to support my explorations.

I use GHC for Haskell, which is the top pick for my language design experiments. I use SBCL and SLIME for Common Lisp hacking, which I don’t do very often these days and only for small exploratory programs. I use MIT Scheme as my Scheme implementation only because that’s what my friend Alexey uses for his language design experiments.

I used to write some Ruby, but I haven’t done any Ruby hacking for quite a while now. I do have some Ruby scripts lying around that I use occasionally (for downloading BibTeX items from MathSciNet and ACM).

Soon I’m going to add Erlang to the list of language I use/play with. One (or both) of Standard ML and Caml are also interesting, primarily because being strict and non-pure makes programming in these language sufficiently different from programming in Haskell.

Well, that’s it folks. There are, of course, other pieces of software, but I don’t use them nearly as often as the above. I have a large collection of music, which I used to listen to using rhythmbox. These days, I prefer to tune in to Radio Paradise. I don’t do any photo editing except for extremely rare case when I have to resize or crop a picture. I use totem for watching videos, but most of the videos I watch are online. I’ve learned to love Emacs dired mode, and don’t start Nautilus, the default file manager in Gnome, except by accident. I don’t play games. I rarely run OpenOffice, mainly when I’m forced to edit/view Microsoft Office files. For quick document preparation, I prefer Google Docs.

Oh, there is a small Windows app that I run under Wine because I haven’t found a good Linux alternative: AP Guitar Tuner.

Non-Confluence of the Perturbative λ-Calculus

Over the last month, I have been trying to prove that the reduction relation defined in my MFPS paper “A Simply Typed Lambda-Calculus of Forward Automatic Differentiation” is confluent. Alas, this is not true. A counterexample is provided by the term \mathsf{T}\,(\lambda x.\;(\lambda y.\;f x y) x), where f is a free variable. There are two reduction sequences that produce two different normal forms:

\displaystyle \begin{array}{cl}   & \lambda x.\;\iota_1\,((\pi_1\,((\mathsf{T}\,f) x)) (\pi_2\,x)) \\ {} + & \lambda x.\;\iota_1\,(\pi_1\,((\mathsf{T}\,(f (\pi_2\,x))) x)) \\ {} + & \lambda x.\;\iota_2 ((f (\pi_2\,x)) (\pi_2\,x)), \\ \neq  & \lambda x.\;\iota_1\,((\pi_1\,((\mathsf{T}\,f) x)) (\pi_2\,x)) \\ {} + & \lambda x.\;\iota_1\,(\pi_1\,((\mathsf{T}\,(f (\pi_2\,x))) (\iota_2\,(\pi_2\,x)))) \\ {} + & \lambda x.\;\iota_1\,(\pi_1\,((\mathsf{T}\,(f (\pi_2\,x))) x)) \\ {} + & \lambda x.\;\iota_1\,((\pi_1\,((\mathsf{T}\,f) (\iota_2\,(\pi_2\,x)))) (\pi_2\,x)) \\ {} + & \lambda x.\;\iota_2\,((f (\pi_2\,x)) (\pi_2\,x)). \end{array}

The latter contains terms of the form \pi_1\,((\mathsf{T}\,M) (\iota_2\,N)), which are semantically zero but don’t reduce. The reduction sequences have been discovered with the help of Redex, which proved to be an extremely valuable tool. An encoding of the perturbative λ-calculus in Redex is available here.

Follow

Get every new post delivered to your Inbox.