In some recent work where I’m building a linear algebraic verison of solution density functional theory, I realized that the algebra was quickly getting out of hand. The work neatly shows how to take traditional mathematical methods for working these problems and turn them into numerical methods.

There’s one problem though. The reviewers didn’t see how its result was different from those ancient mathematical ways. To put their fears to rest, I set out to make a few non-traditional numerical implementations of the theory. What I found was that every time I translate the equations to matrices and vectors, I have to consider lots of special cases so that the code quickly becomes littered with ideas that don’t directly translate to equations in the manuscript. This is bad practice, since you don’t want lots of computations without corresponding, documented definitions in your code. There are 2 ways to solve that problem. Either I could make the manuscript significantly larger, or I could write a compiler taking high-level equations and automagically cranking out implementation code.

The second approach is more elegant, and is exemplified by the FEniCS project. There, code like dE = 0.5*eps*inner(nabla_grad(phi), nabla_grad(phi)) creates a symbolic expression standing in for the electrostatic energy density of a potential field, phi. I say ‘standing in’ for, since the new object, dE just created is a mathematical expression tree – not a real number. The expression tree only says how to do the work. If you want to get a final answer, you have to provide eps and phi on a grid (as a vector or matrix). Then the computer can follow the instructions in the tree step by step to get to the answer.

To see how this can work, I built a simple polynomial class that describes how to multiply and add polynomials of the form:

$$ f(x,y) = \sum_{n,m} c_{n,m} x^n y^m $$

This is a standard form, and can be expanded to a really useful symbolic form by making the definitions:

  • Monomial: \( x^n y^m \)
  • Coefficient: \( c_{n,m} \)
  • Polynomial: Sum of monomials times coefficients.

Each of the definitions above gets represented by its own data structure. Without further ado, my implementation is below:

def applyIter(m, k): # cf Haskell's list monad
    for e in m:
        for j in k(e):
            yield j

# dictionary mapping monomials to coefficient ring
# addition must be implemented for each coefficient
# if there are redundant monomials.
class Polynomial(dict):
    def __init__(self, lst): # merge all keys by addition
        d = {}
        for k, v in lst:
            if d.has_key(k):
                d[k] += v
            else:
                d[k] = v
        dict.__init__(self, d)

    def __mul__(a, b):
        return Polynomial(applyIter(a.iteritems(), \
                    lambda (ai,aj): [(ai*bi, aj*bj) for bi,bj in b.iteritems()]))
    def __add__(a, b): # simple concat
        return Polynomial(applyIter([a.iteritems(), b.iteritems()], \
                    lambda u: u))

class Monomial(tuple):
    def __mul__(a, b):
        return Monomial([i+j for i,j in zip(a, b)])
    def __repr__(self):
        return "Monomial%s"%(tuple.__repr__(self))

def test():
    # 3 x - y
    mx = Monomial([1, 0]) # x^1 y^0
    my = Monomial([0, 1]) # x^0 y^1
    a = Polynomial([(mx, 3.0), (my, -1.0)])
    b = a*a
    print(b)

Running the test() function shows the result:

{Monomial(2, 0): 9.0, Monomial(0, 2): 1.0, Monomial(1, 1): -6.0}

which is the polynomial representing \( 9 x^2 + y^2 - 6 x y\).

If you know the scipy library, you might be asking how this differs from Numpy’s polynomial module. First, I used terminology and concepts from Macaulay2, a special-purpose symbolic algebra package for these kinds of calculations. Second, and more obviously, the representation here is sparse – so that you don’t have to have a vector or matrix of coefficients.

The link with the Macaulay2 project is intentional. In the actual application I’m developing, the monomials are members of a closed, finite symmetry group. Also, the coefficients are radial basis functions rather than numbers. Symbolic algebra with those objects isn’t possible with vanilla polynomials, but will be in my complete code.