Tuesday, January 22, 2013

Writing inference algorithms made easier with NIFTY

Today, my colleagues and I, lead by the heroic efforts of Marco Selig, released a new software package called NIFTY (Numerical Information Field Theory). It's a Python library that is designed to simplify the development, testing, and application of statistical inference algorithms (among other things). It does this in a couple of ways.

First, one generally derives an algorithm with pen and paper (or my favorite, chalk), arriving at an equation or set of equations that can be used to obtain an estimate of some model parameters or continuous field (as in the case of image reconstruction, for example). In our group, we use the mathematical framework of Information Field Theory to derive such equations from first principles. NIFTY was designed to make the translation from mathematical language into computer code as simple and intuitive as possible.

But the real power of NIFTY is that it allows one to write algorithms without (much) consideration of the underlying geometry of the problem. Typically an algorithm that's written to be applied to a 3D cartesian geometry, say, wouldn't work for data that is defined on a sphere. This is not the case with NIFTY because with it handles the geometry, or "space", separately from the operators and fields used to describing the problem. In NIFTY, operators and fields are associated with some space when they are initialized, and will automatically call the appropriate operations for that space when required. This means that one can simply change the space on which objects are defined without fiddling with any of the rest of the code.

The advantages of this approach are that it allows for rapid testing of code. For example, one can develop an idea quickly in 1D and then when the software is stable, very easily apply it later to a 3D problem. It also means that code can be more easily recycled, since, for example, an algorithm designed to solve a problem on the sphere can be reused to solve a problem in some other geometry often by just changing one or two lines of code (redefining the spaces).

Furthermore, NIFTY also provides a lot of helpful tools that are often needed when writing inference algorithms, like transformations between sets of often used conjugate spaces, matrix probing, vector smoothing operations, etc. In all, it aims to make the life of the algorithm developer easier.

To illustrate how NIFTY works (I hope this will be illustrative, anyway), say you have a data vector, d, that is related to some physical signal vector, s, via a simple linear relationship

d = Rs + n

where n is a noise vector and R is an operator (represented here as a matrix) that encodes the effects of your measurement device. If both the signal and noise are well modeled using Gaussian statistics, then it can be shown that the Wiener filter is the optimal estimator for s. The Wiener Filter can be expressed as

m = (S-1 + RN-1R)-1RN-1d

where m is the estimate for s, and S and N are the covariance matrices of the signal and noise, respectively. To be clear, X-1 implies matrix inversion and † implies transpose conjugation.

In NIFTY, each of these vectors and operators will be assigned to special NIFTY Python objects. The vectors are defined as instances of NIFTY "field" classes, and the operators are defined as instances of NIFTY "operator" classes. Each of these has an associated "space" class, which defines the geometry of the problem. Assuming that all of the necessary objects are initialized already, coding the filter is as simple as:

m = D(R.adjoint_times(N.inverse_times(d))

where D, which is the (S-1 + RN-1R)-1 part of the formula, is an instance of the "operator" class with an inverse_times() method defined as

def inverse_times(x):
    return S.inverse_times(x) + \
        R.adjoint_times(N.inverse_times(R(x)))

To apply D when we have only defined D.inverse_times() we generally use the Conjugate Gradient method. That's about it. Somewhere before this a space object has been defined and when R, S, N, and d are initialized the space object is specified. The space definition can be changed to your hearts content, and the code shown here doesn't change a bit.

I think you can see that this code is a very natural translation of the mathematical expressions shown above and furthermore that the distinction between fields, operators and spaces makes the algorithm much simpler to develop and much more flexible to work with.

NIFTY has been released as open source under the GPLv3 license. We have also written an accompanying paper that has been posted to the arXiv today and submitted to the journal IEEE Transactions on Signal Processing. For more information, or to get NIFTY: