This page lists projects which might be suitable for undergraduate research.

Table of Contents¶

From Symbols to Numbers¶

Often we have a set of equations that need to be implemented numerically. To quickly get code working, while ensuring that no mistakes are made, it is often convenient to start with the symbolic expressions and then generate code that performs the numerical computation. This is especially useful when working with minimizers or solvers that can utilize derivative information.

However, when nested deeply, or highly differentiated, even relatively simple expression can become unweildy in their final form. The goal of this project is to provide an efficient and effective set of tools for getting from symbolic expressions to the robust code required to evaluate them numerically.

Use Cases¶

Lagrangian Dynamics¶

Consider the problem of setting up the Euler-Lagrange equations for a triple pendulum. This generates pretty nasty expressions.

Neutron Stars¶

As a use case, we have a project for modelling the physics of neutron stars. Here the star is described by a nuclear equation-of-state (EoS) defined by the function $\mathcal{E}(n_n, n_p)$ which is the energy-density (energy per unit volume) of homogeneous nuclear matter as a function of neutron $n_n$ and proton $n_p$ densities. This expression on its own is rather large as it results from inverting a system of 5 linear equations, but then must be inserted into a further set of equations defining the so-called compressible liquid drop model (CLDM) for nuclei which approximates the construction of a nuclear lattice in the crust of the neutron star by minimizing the total energy including electromagnetic interaction and surface effects. Inserting $\mathcal{E}(n_n, n_p)$ into these equations and taking derivatives as required to define the minimization conditions produces huge expressions that take hours to even simplify.

We practically solved this problem by treating $\mathcal{E}(n_n, n_p)$ as a numerical black box, and decoupling the form of these expressions from the CLDM model. This is not ideal, however, as it required defining an appropriate interface for calling $\mathcal{E}(n_n, n_p)$ which was additional programming work.

Existing Tools¶

Most of our code is in python, so we would like to start there. I am aware of the following tools which might be useful:

• SymPy: This is the "standard" symbolic manipulation package for python. It has quite a few tools, and works reasonably well. Here are some notes, including serious limitations:

  • Simplification: Once large expressions are generated, the simplification tools are very slow. In particular, the default approach expr.simplify() tries many different strategies, then measures the outcomes. Some of these strategies are extremely expensive, so this can literally take hours. Unfortunately, no record is kept about which strategies worked and which did not, so one cannot easily issolate the best approach. One small community project might be to look at the sympy code and at least keep track of which simplification strategy worked the best so the user can call that directly in the future.
  • Numerical Computation: Once an expression is found, fairly efficient numerical implementations can be produced using sympy.lambdify. Performance characteristics of various numerical approaches are discussed here.

  • Caching: Once work is done simplifying and generating numerical functions, it would be great to cache these to disk so they can be reused without future expensive computations. Unfortunately, general expressions and lambdified expressions in particular have no good way to be archived. Expressions can often be converted to strings, but these take some time to process when loading, which can be a performance issue. Caching the lamdified functions is thus preferred from a performance standpoint, but runs into many issues:

    • How to serialize sympy lambdified function?
    • Cannot pickle lambdas from SymPy
    • Python - write symbolic expression (sympy) to txt file
    • Cannot pickle sympy.UndefinedFunction

    Many of these issues have been partially resolved using appropriate settings of advanced pickling libraries such as dill or cloudpickle, but they demonstrate that this approach can be brittle. Another issue is that pickles can "go sour" meaning that if you upgrade your code, or even the underlying python interpreter, the pickles might not be able to be loaded.

    My preferred approach is to use human-readable (and writable) representations and have had reasonable success using sympy.printing.lambdarepr.NumPyPrinter to print an importable module with the following structure:

    from __future__ import division
    import math
    import numpy
    import uncertainties.core
    import uncertainties.unumpy
    import sympy
                                         # Comment to keep nice formatting in 
                                         # notebook... (not sure why needed)
    class Mixin(object):
    _names = {names}
    
    @staticmethod
    def {name}(*v):
        MutableDenseMatrix = numpy.array
    
        if any([isinstance(_v, sympy.Expr) for _v in v]):
            np = sympy
            MutableDenseMatrix = sympy.MutableDenseMatrix
        elif any([isinstance(_v, uncertainties.core.AffineScalarFunc) for _v in v]):
            np = uncertainties.unumpy
        elif any([isinstance(_v, numpy.ndarray) for _v in v]):
            np = numpy
        else:
            np = math
    
        exp, log, pi = np.exp, np.log, numpy.pi
        {args} = v
        return {expr}
    

    This provides a human friendly form of importable code with minimal performance hits. It also allows the functions to be called symbolically, with numpy, or with numbers. However, it is a bit of a hack and requires careful attention to where these "cache" files are stored.

• Theano: This project is designed to produce optimized code for symbolic expressions using the CPU, multiple threads, and GPU accelleration if possible. It works, supports CSE, and is quite fast, but unfortunately, the expression manipulation is all done internally and the results are not easily extracted. It is also a somewhat large package and can be difficult to install (as opposed to sympy which is pure python). I also find Theano to be less intuitive to use than SymPy.

  In 2013 people looked at translation between SymPy and Theano. Here are some discussions about using both SymPy and Theano. This does not seem to be currently an active topic:

  • Theano + Sympy for system of ODEs
  • theano_sympy github project.
  • Matthew Rocklin's blog] about Sympy and Theano (2013)
    • Code Generation
    • Scalar Simplification
    • Matrix Expressions
    • Using SymPy within Theano

The crux of this project thus seems to be about combining an efficient symbolic computation path for specifying non-linear equations with a CSE aware simplification process for computing derivatives etc. and finally generating somewhat optimized numerical code. (Using Numba would be great too.)

Numerical Differentiation¶

In many applications, one needs to compute derivatives numerically. In principle this is a straightforward application of the formula:

$$ f'(x) = \frac{f(x + h) - f(x - h)}{2h} + h^2\frac{f'''(x)}{6} + \order(h^4). $$

However, in practice, one must trade truncation errors (minimized by making $h$ as small as possible) with roundoff errors amplified in the subtraction. To be explicit, assume that $f(x)$ can be computed to a relative precision of $\epsilon$. The roundoff error in this expression is thus of the order:

$$ \text{roundoff error (relative)} \sim \frac{\epsilon}{\sqrt{2}h}. $$

Balancing these we obtain an optional value of:

$$ h_c \approx \sqrt[3]{\frac{6\epsilon f(x)}{\sqrt{2}f'''(x)}} \approx 1.6x\sqrt[3]{\epsilon} $$

assuming that $f'''(x)/f(x) \sim x^{-3}$. Even with this optimal step size, the relative error will be of order $\epsilon^{2/3}/2x$ which for double precision reduces the precision from 16 digits to 10. Here we demonstrate these fomulae these for the function $\sin'(x) = \cos(x)$ at $x=1$:

For many functions, the accuracy $\epsilon$ is far less than machine precision (for example, $f(x)$ might involve running a solver where a low tolerance is required for reasonable performance). A further complication is that the behavior of the function may not be known in advance, making the choice of optimal $h$ tricky because of the factor. Thus, although fast, in general this approach is not ideal.

Richardson Extrapolation¶

One might wonder if higher accuracy can be obtained when the function can be computed to high (machine) precision. The answer is yes and lies in applying the technique of Richardson Extrapolation - sometimes called "The deferred approach to the limit". The idea is to compute the finite difference formalism for a series of $h > h_c$ where the function is smooth, and then fit the error with a low-order polynomial, extrapolating the result to $h=0$ before roundoff errors spoil the results. Here we use an implementation of this to demonstrate that we can achieve fairly high relative precision (almost machine precision) and compare the actual error with the estimate of the error made while computing the Richardson Extrapolation:

However, this approach fails spectacularly when the function evaluation is noisy.

An interesting point is that the function needs to be computed to about two orders of magnitude more accuracy than the desired accuracy of the derivative (see Roy:2005 and Eca:2009).