lifting bits out of numerical methods

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.