Cleaning Images with Gaussian Process Regression

The choice of covariance function determines the smoothness and the extent of the kernel. In principle, a wide variety of covariance functions could be used, but we will focus on the squared exponential covariance function because of its simplicity and its good approximation to a typical instrumental PSF:

The parameter $a$ in Equation (4) describes the magnitude of the correlation between pixels; it should be compared to the data variance $\sigma^{2}$. The kernel does not depend on $a$ and $\sigma$ individually but on the correlation-to-noise ratio $\alpha\equiv a/\sigma$, due to the fact that $K\propto a^{2}=\alpha^{2}\sigma^{2}$ and $C_{\rm data}\propto\sigma^{2}$, so $(K+C_{\rm data})^{-1}\propto\sigma^{-2}$, canceling $\sigma$ in Equation (3). Hence, in the rest of the paper we will set $\sigma=1$ (so that $C_{\rm data}$ is the identity matrix) and use $a$ to mean $a/\sigma$. The noiseless limit $a/\sigma\rightarrow\infty$ implies a constant, dark image, and a kernel independent of $a$. The other limit, $a/\sigma\rightarrow 0$, implies a noisy image where no useful information about a pixel’s value can be extracted from its neighbors. The parameter $h$ gives the characteristic scale over which the image varies. Larger values of $h$ correspond to broader kernels. The values of the hyperparameters $a$ and $h$, together with a particular choice of the covariance function, parameterize a family of GPR kernels.

Figure 2 shows four examples of GPR kernels with different hyperparameters $a/\sigma$ and $h$. As $a/\sigma$ increases, GPR makes use of more pixels in the region and assigns larger weights to the adjacent pixels, allowing significant negative weights to be used. As $h$ increases, the absolute values of the weights decay more slowly with increasing distance; the weights change over a larger scale. All GPR kernels with a single bad pixel have ring-like structure because the value of the squared exponential covariance function depends only on the Euclidean distance to the central pixel, and they are zero at the center where the bad pixel is located.

The fact that each particular choice of $a$, $h$, and covariance function uniquely determines a convolution kernel has several advantages. One of them is that when there is a continuous region of bad pixels to be assigned zero weights, GPR is able to recompute the weights on all other pixels instead of simply re-normalizing them. To give zero weights to additional bad pixels near the bad pixel to be corrected, all we need to do is remove the corresponding entries in the covariance function and GPR automatically compensates for the loss of information by increasing or reducing the weights of their neighbours. Figure 3 panels (a) through (c) show how additional bad pixels near the central bad pixel re-weight the entire GPR kernel and break its symmetry.

Another advantage of GPR is its straightforward handling of bad pixels on the edges of an image. Similar to the accommodation of additional bad pixels, we only need to remove all points that are outside of the image from the covariance function, and GPR will recompute the weights of the rest of the points. Figure 3 panel (d) gives an example of a truncated GPR kernel due to a bad pixel on the top left corner of the image. Because the bad pixel is at $(1,0)$ relative to the corner, with a kernel width of $9$, we are only able to construct a $5\times 6$ kernel. GPR enables the weights to change accordingly to deal with the loss of information.

To choose the optimal hyperparameters $a$ and $h$, our strategy is to identify a set of well-measured pixels in the image to serve as the training set. Given $a$ and $h$, we can then convolve the entire image by the kernel defined by $a$ and $h$ and compute the subsequent change in value at each pixel in the training set. We find the $a$ and $h$ that minimize the mean absolute deviation of these residuals.

The training set has to be large enough to be representative of the image, and it should focus on the pixels of scientific importance. Generally this means pixels with count rates appreciably higher than the background level. A kernel trained on low-SNR background pixels would simply learn to average a large set of surrounding pixels rather than identifying the correlation between pixels containing astronomically useful flux. A kernel optimized on a high-SNR training set would be good at correcting pixels of scientific importance even for low-SNR images (see Section 5.4). For a typical application, then, the training set should consist of pixels well above the background level but below saturation. Section 4.1 describes astrofix’s choice of training set.

In principle, any covariance function can be used to generate GPR kernels, but only covariance functions that resemble the instrumental PSF would be useful. The squared exponential covariance function is therefore a reasonable choice. One notable alternative that we investigate here is the Ornstein-Uhlenbeck covariance function:

Compared to the squared exponential covariance function, the Ornstein-Uhelenbeck covariance function is not quadratic but linear in Euclidean distance in the exponent. As a consequence, it tends to under-emphasize pixels that are the closest to the central bad pixel, which often contain the most important information. Figure 4 compares the optimized kernels for the star image in Figure 1 using the Squared Exponential covariance function (Equation (4)) with those using the Ornstein-Uhelenbeck covariance function (Equation (5)). The magnitude of the weights in the central region is significantly smaller for the $K_{\rm OU}$ kernels than for the $K_{\rm SE}$ kernels.

One other feature of a $K_{\rm OU}$ kernel is that its more slowly changing function values lead to smaller fluctuations in the weights and less prominent negative weights. As a consequence, most outer pixels have almost zero weights. In fact, for most values of $h$ smaller than the kernel width, almost all of the weight in a $K_{\rm OU}$ kernel is concentrated in the central $5\times 5$ region. Hence, a $K_{\rm OU}$ kernel of width $>5$ is almost identical to a $K_{\rm OU}$ kernel of width $5$, which is not the case for $K_{\rm SE}$ kernels.

As a result of the two aforementioned features, a $K_{\rm OU}$ kernel is usually not as accurate and as flexible as a $K_{\rm SE}$ kernel. Figure 5 shows the histogram of residuals (normalized by the $\rm shot\ \rm noise=\sqrt{\rm count/\rm gain}$) from convolving the star image with the two types of kernels. In terms of the mean absolute deviations of the residuals, the $K_{\rm SE}$ kernel with width 9 performs better than the other kernels, while the $K_{\rm OU}$ kernel with width 9 has almost the same residual distribution as the $K_{\rm OU}$ kernel with width 5.

astrofix does have other covariance functions implemented, but we have found the squared exponential to produce the narrowest residual distribution. We restrict further analysis in this paper to the squared exponential covariance function.

So far, we have discussed GPR kernels that have weights depending only on the Euclidean distance $r$, but for anisotropic images, it is sometimes preferred to use a kernel that has directional dependence. Such a kernel can be particularly useful to the imputation of 2D spectroscopic data. For this purpose, we define the stretched squared exponential covariance function:

which produces kernels determined by three parameters: $a$, $h_{x}$, and $h_{y}$. We assume here that the preferred direction (e.g. the dispersion direction) approximately follows one of the orthogonal axes of the detector, though it would be straightforward to adopt a different angle between the preferred direction of the image and the detector axes. For isotropic images, we expect the optimization to give $h_{x}\approx h_{y}$. For vertical spectroscopic lines (horizontal dispersion) in a slit spectrograph, we expect that $h_{x}<h_{y}$ because the pixel values vary on a larger scale in the $y$-direction (along the slit). Figure 6 compares the isotropic and the stretched $K_{SE}$ kernel optimized for the spectroscopic arc lamp image in Section 5.3, and Figure 10 shows that the stretched kernel produces smaller residuals in this example.

We encourage the user to exercise care with a stretched covariance function. The extra free parameter increases the risk of overfitting the hyperparameters. A stretched kernel can also smooth out a point source on the slit, especially if the hyperparameters are optimized on the sky background, a galaxy stretching the length of the slit, or a lamp flat. Our training image in this case was an arc lamp for wavelength calibration; the best-fit hyperparameters would differ for a galaxy or a twilight spectrum.

The set of pixels used to optimize the GPR should include as many pixels of scientific importance as possible, while avoiding pixels that are close to the background level, saturated, or themselves bad. By default, we include all good (unflagged) pixels that are below $1/5$ of the brightest pixel’s value and $10$ median absolute deviations (MAD) above the median of the background. We take the median and the MAD of the image in lieu of the mean and standard deviation, because of the robustness and convenience of the median. The MAD is an underestimate of the standard deviation: for a Gaussian background distribution, $10$ MAD correspond to about $7\sigma$.

In many cases, the users can replace the default training set selection parameters ($1/5$ and $10$) with better ones by deriving the background distribution and the saturation level from the histogram of the image, but for high-SNR images with sufficient bright pixels, the quality of the training is largely insensitive to the parameters. The default training set usually works sufficiently well without requiring a highly accurate model of the background. For low-SNR images, one has to tune the parameters more carefully or use a high-SNR image taken by the same instrument (and under the same configuration) as the training set.

The training set might include pixels that are themselves good but are close to bad pixels that can ruin the convolution. We remove such pixels during the training process described in the next subsection.

In the high-SNR limit, the training process can be considered as modeling the structure that an instrument imparts to an image, such as the PSF. Therefore, an efficient strategy is to run the training only once for each instrument using a combined training set built from multiple images. However, the optimal kernel may still vary from image to image due to, e.g., variations in the seeing. It is not computationally expensive to redo the training as appropriate for each image (see Appendix A).

In the training process, we convolve the training set pixels with GPR kernels and vary the hyperparameters $a$ and $h$ (or $a$, $h_{x}$, and $h_{y}$) to minimize the mean absolute residual in the training set. We only need to compute the kernel once at fixed values for $a$ and $h$. The computational cost then scales as $N_{\rm train}\times w^{2}$, where $w$ is the width (in pixels) of the kernel and $N_{\rm train}$ is the number of pixels in the training set. We construct an image array of size $N_{\rm train}\times w^{2}$ to perform the convolution.

A problem naturally occurs in the limit of noiseless data, where the optimal value of the parameter $a$ (i.e. $a/\sigma$ in our notation) would be infinite. If $a\rightarrow\infty$ and $h$ is relatively large, evaluating Equation (3) requires the inversion of an ill-conditioned (though not formally singular) matrix. This leads to numerical instability in the constructed GPR kernel, and the optimization procedure can end up in a spurious minimum.

We restrict $a$ to small values by multiplying the mean absolute residual with a penalty function:

Here, $\mu$ can be interpreted as a ”soft” upper bound on $a$. For $a\ll\mu$, $P\approx 1$; for $a>\mu$, $P$ increases exponentially. $\tau$ determines how much $a$ can penetrate the soft upper bound. The best values of $\mu$ and $\tau$ are determined by the largest condition number that the machine can handle in a matrix inversion without much loss of accuracy. Figure 7 shows how the condition number depends on the hyperparameters $a$ and $h$ for 3 different kernel widths. A larger kernel width requires inverting a larger matrix, which leads to a larger condition number. The condition number increases quickly with $h$ for small $h$ and increases quickly with $a$ for larger $h$. The contours indicate the upper bound on $a$ given a desired upper bound on the condition number. For instance, if the maximum acceptable condition number is $10^{8}$, then in the worst case scenario (largest $w$ and $h$), $\mu=2400$ is roughly the upper bound on $a$. By default, we set $\mu=3000$ and $\tau=200$, which corresponds to a maximum condition number of about $6\times 10^{8}$ for $w=9$ and $h=9$. The user can change the upper bound parameters according to their own needs and referring to Figure 7.

Then, for a given image and a chosen training set, astrofix’s steps to tune the GPR hyperparameters are:

1. 1.

   Flag all bad pixels as NaN.

2. 2.

   Construct the image array for pixels in the training set, and remove rows that contain NaN values.

3. 3.

   Propose the values of $a$ and $h$ (or $a$, $h_{x}$, and $h_{y}$), and construct the corresponding kernel.

4. 4.

   Convolve the image with the kernel, restricting the convolution to pixels in the training set.

5. 5.

   Find the mean absolute residual of the convolved image and the original image for pixels in the training set.

6. 6.

   Minimize the mean absolute residual times the penalty function with respect to the parameters.

7. 7.

   Return the optimal parameter values and the residual distribution.

Steps 3 through 5 are wrapped into one single function to be passed to scipy.optimize.minimize, which is responsible for step 6. In step 3, we do not reweight the kernel in the presence of additional bad pixels because of its expense. We simply remove, in step 2, all affected pixels from the training set by checking each row of the image array for NaN values. In step 6, we require that $a\geq 1$ because $a<1$ would imply that the correlation between pixels is smaller than the measurement noise. If $a$ is small, neighboring pixels would contribute little information about a missing data point and our effort would be meaningless. We also require that $1/2\leq h\leq w$. A characteristic scale larger than the kernel width would imply a spatially uniform image, while one much smaller than a pixel would again mean that a pixel’s neighbors do not contain meaningful information to interpolate its value.

The computational expense of the optimization process depends on the size of the training set, but it usually takes only a few seconds. For example, it took $\approx$6 seconds to train on the CHARIS spectrograph in Section 5.2 with $\approx$439,000 pixels in the training set. Appendix A includes a detailed analysis of the computational cost of astrofix.

The optimal kernel found in the training process relies on the assumption that for every bad pixel, no other bad pixels or image edges lie within the $w\times w$ region of the kernel. In this case, correcting the image is equivalent to a simple convolution with the optimal kernel at the location of each bad pixel. However, if there are additional bad pixels or pixels off the edges of the image, we must set their weights to zero and recompute the entire kernel by removing the corresponding entries in the covariance function (as discussed in Section 2.2). When the bad pixel density is large, almost every bad pixel requires computing a new kernel, so the resulting computational cost scales linearly with the bad pixel density and can be expensive (see Appendix A).

The optimal kernel that corresponds to an isolated bad pixel is usually used many times, so we pre-construct it before looping through all bad pixels. Because we correct only one bad pixel at a time and only a small fraction of the image overall, we use the weighted sum of the neighboring pixels rather than convolving the entire image. Finally, we replace the bad pixels with their imputed values only after exiting the loop. This avoids using imputed data as if they were actual measurements.

To sum up, astrofix’s steps to impute the missing data in an image are as follows:

1. 1.

   Pre-construct the optimal kernel for an isolated bad pixel.

2. 2.

   For each bad pixel:

   Check if there are additional bad pixels or image edges nearby.

   If not, use the optimal kernel. Otherwise, construct a kernel that has the largest possible size ($\leq w$ in each direction) and has weights reassigned according to the locations of additional missing pixels.

   Compute the sum of the neighboring pixels weighted by the kernel to find the bad pixel’s imputed value. Store this value.

3. 3.

   Replace the bad pixels with the imputed values.

4. 4.

   Return the fixed image.

Images of star clusters often contain a high density of useful information. While approaches like PSF photometry can handle missing data appropriately (if a good model of the PSF is available), such an analysis is not always possible. Aperture photometry, for example, cannot tolerate missing data. We therefore test astrofix on a real CCD exposure of the globular cluster NGC 104, taken by the SBIG 6303 0.4m telescope at the Las Cumbres Observatory (LCO). We randomly generated 62069 ($\sim$1%) artificial bad pixels and imputed their values.

The left panel of Figure 8 compares the residual of image imputation $\Delta\mu$ by GPR to that by two standard approaches: astropy.convolution (with a Gaussian kernel of standard deviation 1) and the $5\times 5$ median filter used by LACosmic. The residual is normalized by the shot noise. Because most bad pixels that we have randomly selected were near the background level, we only include in the histogram the 2642 bad pixels that originally had counts 10$\sigma$ above the background. This tests astrofix’s ability to restore useful information. The right panels of Figure 8 show the unnormalized residuals in a $25\times 25$ region at the center of the cluster. GPR outperforms the other approaches, especially in dealing with continuous regions of bad pixels. In terms of the outlier-resistant mean absolute deviation, GPR outperforms astropy’s interpolation by about a factor of 2, and median replacement by more than a factor of 3.

Infrared detectors are more specialized than CCDs, and often have a higher incidence of bad pixels. The CHARIS integral-field spectrograph (IFS), for example, uses a Hawaii2-RG array (Peters et al., 2012), and has a bad pixel fraction of about 0.5% (Brandt et al., 2017). Similar instruments include the GPI (Macintosh et al., 2014) and SPHERE (Beuzit et al., 2019) IFSs; these have comparable bad pixel rates. The data processing pipelines for GPI (Perrin et al., 2014) and SPHERE both use box extraction to reconstruct the spectra of each lenslet. This approach cannot accommodate missing data, so both pipelines interpolate to fill in bad pixel values.

We test astrofix on a raw CHARIS image with 26338 intrinsic bad pixels flagged by values of zero. There were 438996 pixels in the training set, and the total computational time was 14.96 seconds. The top panels of Figure 9 show a slice of the raw image as well as the corrected images by GPR, by $5\times 5$ median replacement, and by astropy.convolution ($\sigma=1$). The bottom panels show a slice through the corresponding extracted data cubes. The original data cube comes from the standard CHARIS data reduction pipeline (Brandt et al., 2017) which ignores bad pixels. GPR clearly outperforms both astropy and the median filter and performs comparably well to the standard pipeline. When the data processing pipeline cannot ignore bad pixels (e.g. box extraction for GPI and SPHERE), astrofix would become a useful interpolation method.

The standard GPR approach does not require any special treatment to adapt well to spectroscopic data, but we can apply a stretched GPR kernel to further improve its performance. In a real spectroscopic CCD exposure taken by the FLOYDS spectrograph (Sand et al., 2011) at LCO, we randomly generated 53537 ($\sim$5%) artificial bad pixels, of which 1524 pixels had counts $10\sigma$ above the background. We applied both GPR and stretched GPR (sGPR) to impute the missing data, and their residuals (normalized by the shot noise) are compared in the left panel of Figure 10 with the residuals from astropy.convolution ($\sigma=1$) and the $5\times 5$ median filter. The right panels of Figure 10 show the unnormalized residuals in a $35\times 25$ region. Both GPR and sGPR outperform other approaches by a factor $>4$, and sGPR offers a small improvement over GPR.

The data in this case are from an arc lamp and show little structure perpendicular to the dispersion direction. The results would differ for spectroscopy of a point source. Depending on the strength of spectral features, the image in the dispersion direction might appear smoother than in the slit direction, the reverse of this case. For this reason, we urge caution in the use of sGPR with astrofix.

The previous examples demonstrate astrofix’s ability to repair high-SNR images. In the low-SNR limit, however, the residual minimization process of astrofix favors a uniform kernel with $a\rightarrow 0$ to smooth out the background noise. This reduces astrofix to a simple average, which is suboptimal for restoring the useful information in a pixel. Therefore, we recommend that the user restrict the training to the brighter pixels in the low-SNR image or train on a high-SNR image taken by the same instrument (and under the same configuration) and apply the optimized kernel to the low-SNR image. Because such a kernel is trained on strongly correlated neighbouring pixels, it would be non-uniform and could produce a wider residual distribution for low-SNR bad pixels. However, the optimized kernel would perform substantially better than the uniform kernel in repairing pixels of practical use, such as pixels containing the PSF of a star.

As an example, we compare astrofix and astropy.convolution in terms of PSF photometry and astrometry on the truncated stellar image in Figure 1. We generate artificial bad pixels in its PSF and add white Gaussian noise to lower its SNR. We then apply both astrofix and astropy.convolution to repair the PSF at $\rm SNR=5,6,7,...,40$. We train the GPR kernel on the rest of the image with high SNR and get $a=3.02,h=0.72$. For astropy.convolution, we apply a Gaussian kernel of $\sigma=0.95$ to match the $\sigma$ of the fitted PSF. Finally, we perform PSF photometry and astrometry on the fixed images and compute the deviations of the fitted centroids from the true centroid position ($\Delta x$ and $\Delta y$, in the unit of pixels) and the fractional difference between the fitted flux and the true flux ($\Delta F/F$). For both photometry and astrometry, we assume Gaussian PSFs and use IterativelySubtractedPSFPhotometry and IRAFStarFinder in the Python package Photutils with the default parameters. We repeat the procedure for 1000 random noise samples at each SNR level.

Figure 11 shows the results for two different sets of bad pixels. astrofix produces much smaller systematic errors on all three fitted parameters, showing that the optimized GPR kernel is the more accurate model of the PSF, although, at very small SNRs, it generally produces a slightly larger spread in the fitted parameters. In other words, we trade precision for accuracy by training the kernel on high-SNR pixels, putting greater weights on modeling the PSF instead of averaging out the noise. This strategy maintains the scientific value of astrofix even for low-SNR images, and a high-SNR training set is almost always available from the brighter regions on the same image or other images taken by the same instrument.

A fully uniform background would remove the bias seen in Figure 11. In this case, each neighboring pixel would contain the same information as the bad pixel, and the optimal approach would be a simple average of neighboring good pixels. Our approach with astrofix, however, will better retain small-scale structure atop a smoothly varying background.