Golomb Rulers/Squares/Rectangles

Yesterday I came across an interesting example of the isolation of academic fields from one another.

A common design parameter for antenna arrays is to try to obtain uniform coverage of the aperture plane to get as many independent measurements as possible. Interferometers sample the aperture plane at locations that correspond to the difference vectors between antenna elements. For example, in one dimension, if I put four antennas at locations (0,1,4,6), then that array would sample the difference set of those positions: (1,2,3,4,5,6). However, if I put antennas at (0,1,2,3), the difference set would only include (1,2,3) with 1 occuring 3 times and 2 occuring twice. For the purpose of uniformly sampling an aperture, redundant spacings are lost measurements. We're looking for a minimum-redundancy array.
There have been a few papers in radio astronomy on minimum-redundancy arrays for the more useful 2-dimensional case, including Golay (1971), Klemperer (1974), and Cornwell (1988).

Thinking that this might be a mathematical problem of interest, I ran it by a good friend of mine: Phil Matchett Wood--a mathematician at Princeton. He quickly uncovered the equivalent problem as formulated in math literature: Golomb rulers in 1-D and Golomb rectangles in 2-D. Some relevant papers on the subject are Shearer (1995), Meyer & Jaumard (2005), Robinson (1985), and Robinson (1997). These papers are on the exact problem and describe applications to "radar and sonar signal detection", obviously referring to the need for independent aperture samples in radar and sonar interferometers. Somehow, the differing nomenclature between these fields was never quite bridged, and so there has not been and cross-referencing between these two formulations of the same underlying problem. People like to talk about the possibility that relevant research in one field goes unnoticed by other fields. This is the first time I've across it myself, though.

More on MCMC in Python

Markov-Chain Monte Carlo (MCMC) seems to be a promising technique for the calibration/imaging problem that we are facing with our experiment the Precision Array for Probing the Epoch of Reionization (PAPER). Yesterday, in addition to taking a crash-course in MCMC, I also started playing with PyMC, which implements, among other things, MCMC using Metropolis-Hastings chains. A first shot at a simple fitter using PyMC went something like this:

import pymc, numpy as n, pylab as p
from pymc import deterministic

x = n.arange(-10., 10, .01)

def actual_func(a, b, c): return a*x**2 + b*x + c

sig = .1
tau = 1/sig**2
noise = n.random.normal(0, scale=sig, size=x.size)
mdata = actual_func(1., 2., 3.)
mdata += noise

a_fit = pymc.Uniform('a', 0., 2.)
b_fit = pymc.Uniform('b', 1., 3.)
c_fit = pymc.Uniform('c', 2., 4.)
tau_fit = pymc.Uniform('tau', tau/3, 3*tau)

@deterministic
def func(a=a_fit, b=b_fit, c=c_fit): return actual_func(a, b, c)

mdata_var = pymc.Normal('mdata', mu=func, tau=tau_fit,
    value=mdata, observed=True)

mdl = pymc.Model([a_fit, b_fit, c_fit, d_fit, mdata_var, tau_fit])
mc = pymc.MCMC(mdl)
mc.sample(iter=1000,burn=250)  

a_hist, a_edges = n.histogram(a_fit.trace(), 40, (0,2))
a_bins = (a_edges[:-1] + a_edges[1:])/2
p.plot(a_bins, a_hist)

p.show()

I chose this example over the ones that came with PyMC because it was much closer to the kind of calibration problems we will be trying to solve with PAPER. An interesting point that came up early on was the reliance of MCMC on an estimate of noise levels in the data. I remember this from the Maximum Entropy Method (MEM) for deconvolution that I coded up in AIPY. An interesting technique, used here, is to actually include the noise level (tau_fit) as a variable that is fit for. This way you can specify a starting guess for noise, along with a confidence interval, and not have to pull a magic number out of a hat. In this example, the fitter does a good job of accurately determining noise levels. I think what happens in this code is that the current guess for the noise level is used as part of the calculation that determines the next state to jump to in the Markov chain, and that new state make include a revision of the noise level. This clearly might be instable for more complex systems, so I imagine some amount of care must be exercised in leaving noise as a free parameter.

Image Deconvolution with Markov-Chain Monte Carlo

I've decided to morph this blog to be more about short updates pertaining to current ideas I'm thinking about, rather than the long-winded philosophical rants I've posted so far. Hopefully this might keep me more engaged as a blogger and maybe even help me keep better track of the things I'm working on. Starting up this blog again, by the way, is a shameless procrastination technique, since my dissertation is due in about 1 month, and all my writing energy should really be focused on that...

After attending an SKA Imaging Workshop in Socorro, NM a couple of weeks ago, I've developed an interest in Bayesian statistics and Markov-Chain Monte Carlo (MCMC) techniques as they pertain to interferometric imaging. Having never taken a stats course, I'm scrambling a little to absorb the vocabulary I need to understand papers written on the subject. Fortunately, in this era of wikipedia, getting up to speed isn't that hard. After reading wiki articles on MCMC, Markov Chains, and the Metropolis-Hastings algorithm, I dived into EVLA Memo 102, which talks about a first shot at using MCMC for image deconvolution.

The Maximum Entropy Method (MEM) is a classic deconvolution technique (one I've already reimplemented for AIPY), but I'd like to go a bit further down this road. According to the standard implementation (which I gleaned from reading Cornwell & Evans (1984) and Sault (1990)) this algorithm uses a steepest descent minimization technique based on the assumption of a nearly diagonal pixel covariance matrix (i.e. the convolution kernel is approximately a single pixel). While this is an effective computation-saving assumption, I found that for the data I was working with, this assumption lead to the fit diverging when I started imaging at finer resolutions.

I think MCMC, by not taking the steepest decent, might be able to employ the diagonality assumption more robustly. I also think it's high time that deconvolution algorithms make better use of priors. The spatial uniformity prior in MEM makes it powerful for deconvolving extended emission, while the brightest-pixel selection technique in CLEAN makes it effective for deconvolving point sources. There's no reason we can't build a prior explicitly for a deconvolution algorithm that tells it to prefer single strong point sources over many weaker point sources, but also tells it that when all else is equal, entropy should be maximized.