Nonlinear fitting of undersampled discrete datasets in astronomy

Abstract: Data analysis and interpretation often relies on an approximation of an empirical dataset by some analytic functions or models. Actual implementations usually rely on a non-linear multi-dimensional optimization algorithm, typically Levenberg--Marquardt (LM) or other flavors of Newtonian gradient methods. A vast majority of datasets in optical and infrared astronomy are represented by values on a discrete grid because the actual signal is sampled by regularly shaped pixels in the light detectors. Here we come to the main problem of nearly all widely used implementations of nonlinear optimization methods: the function that is being fitted is evaluated at central pixel positions rather than integrated over the pixel areas. Therefore, the best-fitting set of parameters returned by the minimization routine might not be the best representation of the observed dataset, especially if a dataset is undersampled. For example, a central pixel of a 1D Gaussian with a dispersion of 1 pix (2.36 pix FWHM; so not too strongly undersampled) will be about 4.2% lower than its central evaluated value if integrated. To handle this effect properly, one needs to perform numerical or analytic integration of a model within the pixel boundaries. We will discuss possible computationally efficient solutions and test our preliminary implementation of a nonlinear fitting using LM minimization that correctly accounts for the discrete nature of the data.