Multiband Probabilistic Cataloging: A Joint Fitting Approach to Point-source Detection and Deblending

Within the framework of probabilistic cataloging, catalogs are best thought about as samples drawn from the model space—many of the challenges of cataloging are common to inference in high-dimensional spaces. We define a catalog as a sample containing N sources ${\{({x}_{n},{y}_{n},{f}_{n},...)\}}_{n\,=\,1}^{N}$ describing some data. We only consider point sources in this paper, so our catalog will typically look like ${\{({x}_{n},{y}_{n},{f}_{1,n},{f}_{2,n},...,{f}_{B,n})\}}_{n\,=\,1}^{N}$, where ${\{{f}_{b,n}\}}_{b=1}^{B}$ represent the fluxes in B different bands for source n and $({x}_{n},{y}_{n})$ denotes the source position in a fiducial image. We assume the astrometric transformations to (x, y) in other images to be known and fixed. Alternative catalogs that include extended sources such as galaxies require more complicated parameterizations to describe the orientation and morphology of each source.

Adapting notation from Portillo et al. (2017), we express the expected counts at each pixel grid coordinate $(x,y)=(l,m)$ for band b as ${\lambda }_{{lm}}^{b}$, with units of photoelectrons. This is calculated as a sum of the background in band b, ${I}_{\mathrm{sky}}^{b}$, and the contributions of nearby sources:

$$\begin{eqnarray}&&{\lambda }_{{lm}}^{b}={I}_{\mathrm{sky}}^{b}+\sum _{n=1}^{N}{f}_{b,n}{{ \mathcal P }}_{b}(l-{x}_{b,n},m-{y}_{b,n}).\end{eqnarray} \tag{ 1 }$$

In the above equation, ${{ \mathcal P }}_{b}({\rm{\Delta }}x,{\rm{\Delta }}y)$ is the pixel-convolved PSF extracted by the standard SDSS pipeline for the center of our image. For our observations, ${\lambda }_{{lm}}^{b}$ is large enough that the expected noise is approximately Gaussian with a standard deviation of $\sqrt{{\lambda }_{{lm}}^{b}}$ photoelectrons. For data ${k}_{{lm}}^{b}$ over n_b images with width W and height H, the likelihood is

$$\begin{eqnarray}&&{ \mathcal L }=\prod _{b=1}^{B}\prod _{l=1}^{W}\prod _{m=1}^{H}\displaystyle \frac{1}{\sqrt{2\pi {\lambda }_{{lm}}^{b}}}\exp \left(-\displaystyle \frac{{\left({k}_{{lm}}^{b}-{\lambda }_{{lm}}^{b}\right)}^{2}}{2{\lambda }_{{lm}}^{b}}\right).\end{eqnarray} \tag{ 2 }$$

For the purpose of MCMC sampling, we calculate the log-likelihood, which turns our products into sums over pixels and bands:

$$\begin{eqnarray}&&\mathrm{log}{ \mathcal L }=\sum _{b=1}^{B}\sum _{l=1}^{w}\sum _{m=1}^{H}-\displaystyle \frac{1}{2}\mathrm{log}{\lambda }_{{lm}}^{b}+c-\displaystyle \frac{{\left({k}_{{lm}}^{b}-{\lambda }_{{lm}}^{b}\right)}^{2}}{2{\lambda }_{{lm}}^{b}}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\approx \sum _{b=1}^{B}\sum _{l=1}^{w}\sum _{m=1}^{H}-\displaystyle \frac{{\left({k}_{{lm}}^{b}-{\lambda }_{{lm}}^{b}\right)}^{2}}{2{\lambda }_{{lm}}^{b}}.\end{eqnarray} \tag{ 4 }$$

When sampling, we only need to compute delta log-likelihoods ${\rm{\Delta }}\mathrm{log}{ \mathcal L }$, so we are at liberty to discard the additive constant c. While we assume a Gaussian likelihood on the pixels of the images, we choose to also remove the logarithmic term $-\tfrac{1}{2}\mathrm{log}{\lambda }_{{lm}}^{b}$ when computing ${\rm{\Delta }}\mathrm{log}{ \mathcal L }$. Because we use ${\lambda }_{{lm}}^{b}$ as an estimator of the pixel-level variance as shown in Equation (2), an imperfectly estimated PSF will misestimate source fluxes and thus ${\lambda }_{{lm}}^{b}$, meaning that our estimate of the pixel-level variance will be incorrect. Furthermore, the logarithmic term is more costly to compute and has little effect compared to the last term as long as ${\lambda }_{{lm}}^{b}\gg 1$. The range of ${\lambda }_{{lm}}^{b}$ will depend on the background and noise level—in the case where the background level is high enough and the noise is small compared to the background, the above condition should hold.

The Bayesian element of catalog inference allows us to impose a number of priors on the parameters of interest. These priors, used in combination with the log-likelihood during sampling, are outlined in Appendix A.

From the assumption that all point sources in our model are independently realized and neither spatially nor spectrally correlated with each other, one can construct a hierarchical Bayesian model. In a hierarchical model, each point source has priors on its parameters (e.g., flux, color), which are then conditioned by hyperparameters that can also float as free parameters in the fit. These hyperparameters might, for example, parameterize the flux distribution, color distribution, or spatial distribution of a source population. Constraints on these hyperparameters (i.e., hyperpriors) are analogous in function to priors on source parameters and are part of what makes the model hierarchical. When calculating marginalized estimates on catalog source properties, one can also marginalize over hyperparameter uncertainties. Hierarchical models can be abstracted to any level, with hyperparameters of hyperparameters ad infinitum, though in practice the usefulness of models with multiple levels of hyperparameters depends on the quality of data. Nonetheless, the hierarchical Bayesian framework of probabilistic cataloging exploits the ability to propagate prior knowledge from any number of parameters to later parts of the analysis, where those parameters can be constrained with more comprehensive uncertainties. A probabilistic graphical model of our hierarchical model is shown in Figure 1. The hyperparameters in our model are fixed in the following analysis, though in principle they can be floated.

Zoom In Zoom Out Reset image size

One issue that arises when cataloging an image concerns the fact that the number of unknowns (i.e., the number of sources) is an unknown itself. From a modeling perspective, the challenge lies in the fact that each source has a number of associated parameters, so models with varying numbers of sources have varying dimensionality. Furthermore, sampling across models of different dimension (${ \mathcal D }(x)\ne { \mathcal D }(x^{\prime} )$) would not be bijective, meaning that the Markov chain would not satisfy detailed balance. A solution to this issue is provided by Green (1995) through a modified algorithm known as Reversible Jump MCMC (RJ-MCMC). RJ-MCMC samples across models by drawing random auxiliary variables u and $u^{\prime}$ with densities g(u) and $g(u^{\prime} )$ to match the dimensions of the initial and final states, where

$$\begin{eqnarray}&&{ \mathcal D }(x)+{ \mathcal D }(u)={ \mathcal D }(x^{\prime} )+{ \mathcal D }(u^{\prime} )\end{eqnarray} \tag{ 5 }$$

and ${ \mathcal D }$ is the dimension operator (Daylan et al. 2017; Green 1995). In traversing between model spaces of different dimension, $x^{\prime}$ is determined by $x,u$ and x is determined by $x^{\prime} ,u^{\prime}$. So long as the mapping $(x,u)\leftrightarrow (x^{\prime} ,u^{\prime} )$ is a diffeomorphism (i.e., bijective and differentiable in both directions), one can construct proposals so that RJ-MCMC has the same convergence properties as traditional MCMC methods. In the context of probabilistic cataloging, these auxiliary variables are related to model sources being added or removed (often referred to as "births" and "deaths"). For example, in the case of a birth, u contains the position and flux of the new source, and $u^{\prime}$ is the empty set. Sources may also be split or merged together with similar auxiliary variables.

While one can explore the parameter space of a catalog model with fixed dimension, the transdimensional nature of probabilistic cataloging enables us to explore a generalized catalog space, which we define as the union of fixed-dimensional catalog parameter spaces. This is written as

$$\begin{eqnarray}&&{ \mathcal C }=\bigcup _{N={N}_{\min }}^{{N}_{\max }}{{ \mathcal C }}_{N}=\bigcup _{N={N}_{\min }}^{{N}_{\max }}{X}_{N}\times {Y}_{N}\times {{ \mathcal F }}_{N}\times ...,\end{eqnarray} \tag{ 6 }$$

where N is the number of sources in a given catalog model. So long as detailed balance is satisfied, samples drawn from the stationary distribution of the Markov chain inherently capture covariances across models of differing dimensionality (differing numbers of sources) by exploring the posterior distribution of N. This enables direct, statistically consistent inference on the number of observed sources in an image.

Our implementation of probabilistic cataloging uses Metropolis–Hastings sampling to converge on the stationary distribution. We provide the acceptance factors computed for various proposals in Appendix B.

The code used in this work is an extension of Lion, a significantly faster implementation of PCAT. Here we describe a few of the performance improvements in the updated implementation.

The most computationally expensive part of each step in the Markov chain is constructing and evaluating the model image for a proposed catalog. For a 100 × 100 pixel single-band image in a crowded field, the implementation of PCAT used in Portillo et al. (2017) took ∼10⁹ model evaluations over 12 CPU-days to converge. That implementation evaluates the model image by looping over sources and then looping over pixels, using bilinear interpolation on the PSF to determine the contribution of the current source to the current pixel. Lion loops over catalog difference entries instead, saving significant computational time. However, constructing the model image remains the dominant computational cost.

To speed up model evaluation, Lion rewrites PSF subpixel shifts as a matrix multiplication operation, which allows us to take advantage of highly optimized libraries such as CBLAS.⁵ Consider the model image for a point source i at location (x_i, y_i). We would like to compute the model image in a postage stamp around the pixel containing the source, choosing a postage stamp large enough to contain the flux from the source but not so large as to introduce needless computation. For each pixel j in the postage stamp, the pixel-convolved PSF must be evaluated at the center of the pixel $({x}_{j},{y}_{j})$:

$$\begin{eqnarray}&&{\lambda }_{i,j}={f}_{i}{ \mathcal P }({x}_{j}-{x}_{i},{y}_{j}-{y}_{i}).\end{eqnarray} \tag{ 7 }$$

The PSF can be described by functions for each pixel in the postage stamp that only depend on (Δx_i, Δy_i), the separation between the source and the center of the pixel it is in:

$$\begin{eqnarray}&&{\lambda }_{i,j}={f}_{i}{{ \mathcal P }}_{j}({\rm{\Delta }}{x}_{i},{\rm{\Delta }}{y}_{i}).\end{eqnarray} \tag{ 8 }$$

These functions ${{ \mathcal P }}_{j}$ are approximated as third-degree polynomials in Δx and Δy:

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal P }}_{j}({\rm{\Delta }}x,{\rm{\Delta }}y) & \approx & {q}_{j,1}+{q}_{j,2}{\rm{\Delta }}x+{q}_{j,3}{\rm{\Delta }}y+{q}_{j,4}{\rm{\Delta }}{x}^{2}\\ & + & {q}_{j,5}{\rm{\Delta }}x{\rm{\Delta }}y+{q}_{j,6}{\rm{\Delta }}{y}^{2}+{q}_{j,7}{\rm{\Delta }}{x}^{3}\\ & + & {q}_{j,8}{\rm{\Delta }}{x}^{2}{\rm{\Delta }}y+{q}_{j,9}{\rm{\Delta }}x{\rm{\Delta }}{y}^{2}+{q}_{j,10}{\rm{\Delta }}{y}^{3}.\end{array}\end{eqnarray} \tag{ 9 }$$

For a fixed PSF, the coefficients q_j,k can be precomputed. Then, if a matrix ${\boldsymbol{Q}}$ is constructed with ${{\boldsymbol{Q}}}_{{jk}}$ = q_j,k and the column vector

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{X}}}_{i} & = & {f}_{i}\times (1,{\rm{\Delta }}{x}_{i},{\rm{\Delta }}{y}_{i},{\rm{\Delta }}{x}_{i}^{2},{\rm{\Delta }}{x}_{i}{\rm{\Delta }}{y}_{i},{\rm{\Delta }}{y}_{i}^{2},\\ & & \qquad {\rm{\Delta }}{x}_{i}^{3},{\rm{\Delta }}{x}_{i}^{2}{\rm{\Delta }}{y}_{i},{\rm{\Delta }}{x}_{i}{\rm{\Delta }}{y}_{i}^{2},{\rm{\Delta }}{y}_{i}^{3})\end{array}\end{eqnarray} \tag{ 10 }$$

is constructed for the source, the postage stamp for the source (flattened into a column vector) is simply the matrix product ${{\boldsymbol{QX}}}_{i}$. Collecting the column vectors for all sources ${{\boldsymbol{X}}}_{i}$ into a single matrix ${\boldsymbol{X}}$, all of the postage stamps are columns in the matrix product ${\boldsymbol{QX}}$. Once this matrix product is computed, the postage stamps can be extracted and pasted onto the larger model image. This procedure agrees with bilinear interpolation to ∼0.1% and takes less than 1 ms for 500 sources in a 100 × 100 pixel image, compared to 100 ms for the previous implementation.

In addition to speeding up model evaluation, tuning the step sizes and evaluating multiple proposals with one model image makes each step in the chain more efficient so that only 10⁶ steps are required. Given the same single-band image, the updated implementation takes only half an hour on one CPU to converge, over 500 times faster than the previous implementation. Other than its drastically improved run time, this implementation of probabilistic cataloging shares many of the advantages and disadvantages of the previous implementation. More information about the updated implementation, Lion, will be detailed in a future paper. The code used in this work, Lion 0.1, is publicly available on GitHub (https://github.com/RichardFeder/multiband_pcat).

We first generate a mock SDSS data set in three bands (r, i, g). Our point sources have r-band fluxes drawn from a flux distribution with a power-law slope α = −2, after which we draw colors $(r-i)\sim { \mathcal N }(0.25,0.5)$ and $(g-r)\sim { \mathcal N }(0.25,0.5)$ to calculate i- and g-band fluxes. We fix the background levels in our bands to ${\boldsymbol{b}}=({b}_{r},{b}_{i},{b}_{g})=(179,310,115)$ in analog–digital units (ADU). These values are based off of empirical background estimates from SDSS run 2583, camcol 4, field 136, which we will use in Section 6. We fix the PSF in all bands to the SDSS r-band PSF template from the same field. This PSF is used to generate our mock catalog sources and is assumed to be known perfectly in our mock tests. For our first mock runs, we assume the astrometric solution of each source to be perfect across all bands by setting (x, y)_r = (x, y)_i = (x, y)_g. Once the model image is generated in each band, we add Gaussian noise proportional to the square root of the variance, where the variance in each image is obtained in ADU by dividing the model image by a gain set to 4.62 across all bands.

Figure 2 shows the completeness and false discovery rates for the catalog ensembles of one-, two-, and three-band PCAT runs. Because PCAT gives us an ensemble of catalogs, the completeness we report is an average over catalog samples. For each sample, we match a sample source with a true source if the positions of the two sources differ by less than 0.5 pixels and their r-band magnitudes are within 0.5 mag of each other. The choice of these matching criteria is explained in Section 6.2. We are primarily concerned with understanding the improvement that comes from adding multiple bands to the fit.

Zoom In Zoom Out Reset image size

We can estimate the improvement in completeness by calculating how the magnitude threshold for detecting a source changes with additional observations. Given a single isolated source with flux f in a uniform background B, we wish to calculate $f^{\prime}$, the flux that yields the same S/N by combining n observations:

$$\begin{eqnarray}&&{\left({\rm{S}}/{\rm{N}}\right)}^{2}=\displaystyle \frac{{f}^{2}}{{N}_{\mathrm{eff}}B+f}=\displaystyle \frac{{nf}{{\prime} }^{2}}{{N}_{\mathrm{eff}}B+f^{\prime} }.\end{eqnarray} \tag{ 11 }$$

Taking the approximation N_eff B ≫ f, this reduces to

$$\begin{eqnarray}&&f^{\prime} =\displaystyle \frac{f}{\sqrt{n}}.\end{eqnarray} \tag{ 12 }$$

For f = 250 ADU (r ≈ 22.2), we calculate ${\rm{\Delta }}m=m^{\prime} -m$ to be 0.37 and 0.60 for n = 2 and n = 3 observations, respectively. As seen in Figure 2, the improvements in catalog completeness by adding i and g band to runs on mock data are consistent with these predictions. By comparing the two runs on three-band data (green and red lines), one can see that the color prior plays an important role in enhancing point-source detection past r ∼ 21 and reducing false positives.

In practice, astrometric calibration can be compromised by a number of instrumental and other systematics. To understand how astrometric miscalibration affects the multiband catalog ensemble, we take the same 100 × 100 pixel mock SDSS data and perturb what would otherwise be perfect astrometry. We consider a two-band case (r and i), where i-band model source positions are artificially perturbed.

Figure 3 shows the catalog ensemble completeness and false discovery rates of multiband PCAT for cases with offsets $\delta x=0.005,0.01,0.02,0.05$, and 0.1 pixels. Subpixel perturbations of bright sources bias the delta log-likelihood of a sampling proposal. For a Gaussian problem, this can be approximated as ${\rm{\Delta }}\mathrm{log}{ \mathcal L }\sim {(\delta x/2{\sigma }_{x})}^{2}$, where σ_x is the uncertainty in a source's position. In these cases the likelihood may favor a scenario where bright sources (i.e., sources with small σ_x) are split into pairs of overlapping sources. Indeed, the brightest sources in our mock data set get oversplit when the astrometry is perturbed at the level of 10⁻² pixels, and the effect becomes more severe for larger miscalibrations.

Zoom In Zoom Out Reset image size

Because the two sources are generated at small separations, the condensed catalog procedure from Portillo et al. (2017) is not feasible for making associations between catalog samples and true sources. Instead, a modified procedure is used to associate sample sources with the two mock sources. For consistency, we refer to the fainter of two sources as source 1 and the brighter as source 2. For each sample:

• 1.  
  If the sample contains one source:

  • (a)  
    If ${f}_{2}/{f}_{1}\gt 1$, then associate the sample with the brighter source, so long as its position does not deviate from that of the bright source by more than the separation of the two real sources. This is motivated by the fact that samples with a single source are likely blended versions of the two real sources, from which an association is stronger with the brighter source.
  • (b)  
    If ${f}_{2}/{f}_{1}=1$, associate each sample source with the closest mock source.
• 2.  
  If the sample contains two or more sources, determine the source closest to source 2 and make that association if it is within the separation offset s of source 2. Then do the same with the closest remaining sample source for source 1.

In this procedure not all sample sources make associations, either because their positions deviate too far from the real sources or because they are in a catalog sample containing more than two sources and were not associated first. This procedure has not been optimized fully, but it should provide adequately robust source associations for the purposes of this study.

To assess the performance of multiband PCAT in these various scenarios, the n-source prevalence p_n (i.e., the fraction of catalog samples with source number n) is calculated for n = 1, n = 2, and n ≥ 3. For a converged chain, the prevalence may be interpreted as the posterior probability for a catalog with N_star = n. While in Portillo et al. (2017) the term is used to refer to individual sources in a condensed catalog, the n-source prevalence here is label invariant. A similar source prevalence analysis can be found in Jones et al. (2015), though our analysis does not explore the effect of relative background levels on the fit. In all runs, p_>2 was measured to be ≤0.01. As such, two-source prevalences p₂ are shown in Figure 5 as a function of separation offset, from which p₁ ≈ 1 − p₂.⁶

Zoom In Zoom Out Reset image size

Figure 5 compares the two-source prevalence as a function of separation between one-band (r), joint multiband (r, i, g), and co-add runs. As expected, the co-add and joint multiband PCAT runs detect two sources more often than corresponding single-band runs. As a check on PCAT's performance, we estimate the source separation that is needed to deblend two neighboring sources. The one- and two-source prevalences are related by the ratio of their Bayesian evidences:

$$\begin{eqnarray}&&\displaystyle \frac{{p}_{2}}{{p}_{1}}=\displaystyle \frac{{Z}_{2}}{{Z}_{1}}=\displaystyle \frac{\pi (n=2)}{\pi (n=1)}\displaystyle \frac{\int { \mathcal L }({\theta }_{1},{\theta }_{2})\pi ({\theta }_{1},{\theta }_{2})d{\theta }_{1}d{\theta }_{2}}{\int { \mathcal L }({\theta }_{1})\pi ({\theta }_{1})d{\theta }_{1}},\end{eqnarray} \tag{ 13 }$$

where θ_1,2 is the set of parameters describing the first or second star. For this derivation, we assume that the two-source likelihood is strongly peaked around the true parameters and is approximately an independent Gaussian for each parameter, ignoring possible covariance between parameters. We also approximate the prior as flat in the region where the likelihood is peaked, with an effective width equal to the reciprocal of the value of the prior at the true parameters ${\theta }_{1,\,2}^{* }$. Then, the evidence integral can be approximated as

$$\begin{eqnarray}&&\mathrm{ln}{Z}_{2}\approx \mathrm{ln}\pi (n=2)+\mathrm{ln}{ \mathcal L }({\theta }_{1}^{* },{\theta }_{2}^{* })+\sum _{\theta \in \{{\theta }_{1},{\theta }_{2}\}}\mathrm{ln}\displaystyle \frac{{\sigma }_{\theta }}{{w}_{\theta }},\end{eqnarray} \tag{ 14 }$$

where σ_θ is the uncertainty on parameter θ, estimated from the second derivative of the log-likelihood, and w_θ is the effective width of the prior. This approximation of the log evidence can be interpreted as the sum of the log prior for n = 2, the maximum log-likelihood, and the log ratio of the uncertainty over the prior width for each parameter (these log ratios are always negative since the uncertainty must be smaller than the prior width). For the one-source likelihood, we assume that the likelihood is strongly peaked around parameters $\tilde{\theta }$ corresponding to the total flux and center of flux, giving

$$\begin{eqnarray}&&\mathrm{ln}{Z}_{1}\approx \mathrm{ln}\pi (n=1)+\mathrm{ln}{ \mathcal L }(\tilde{\theta })+\sum _{\theta \in {\theta }_{1}}\mathrm{ln}\displaystyle \frac{{\sigma }_{\theta }}{{w}_{\theta }}.\end{eqnarray} \tag{ 15 }$$

We approximate the delta log-likelihood between the true two-source parameters $({\theta }_{1}^{* },{\theta }_{2}^{* })$ and the blended source parameters $\tilde{\theta }$ by calculating

$$\begin{eqnarray}&&\mathrm{ln}{ \mathcal L }({\theta }_{1}^{* },{\theta }_{2}^{* })-\mathrm{ln}{ \mathcal L }(\tilde{\theta })\approx \sum _{b=1}^{{n}_{\mathrm{bands}}}\sum _{l=1}^{w}\sum _{m=1}^{H}\displaystyle \frac{{\left({\rm{\Delta }}{\lambda }_{{lm}}^{b}\right)}^{2}}{2{\lambda }_{{lm}}^{b}},\end{eqnarray} \tag{ 16 }$$

where ${\rm{\Delta }}{\lambda }_{{lm}}^{b}$ is the difference between the two-source and blended source model images and ${\lambda }_{{lm}}^{b}$ is the two-source model image. These approximations neglect the possible covariance between the noise and the maximum likelihood two-source and one-source model parameters. The differences in the prior terms $\mathrm{ln}({\sigma }_{\theta }/{w}_{\theta })$ contribute significantly to the difference in evidence: they sum to ≈22, with the x, y terms contributing the most at ≈6 each. The flux and color terms contribute slightly less (≈4 and 3 each, respectively), while the prior on the number of sources contributes 0.5 for each additional degree of freedom (see Section A.2). Using p₁ ≈ 1 − p₂ and denoting ${\rm{\Delta }}\mathrm{ln}Z=\mathrm{ln}{Z}_{2}-\mathrm{ln}{Z}_{1}$,

$$\begin{eqnarray}&&{p}_{2}\approx \displaystyle \frac{{e}^{{\rm{\Delta }}\mathrm{ln}Z}}{1+{e}^{{\rm{\Delta }}\mathrm{ln}Z}},\end{eqnarray} \tag{ 17 }$$

so we report the separations where Δln Z = 4.6, corresponding to an expected p₂ = 0.99. These threshold source separations (also shown in Figure 5) are consistent with near-unity two-source prevalences obtained using PCAT for both single-band and multiband runs.

The multiband runs done on r-, i-, and g-band data yield higher prevalences at small separations than corresponding runs on co-added data. This is evidence that imposing color priors can improve the fit of the resulting catalog ensemble. Note that the color priors used in these runs are fairly broad and uninformative. The degree of improvement one might expect from spectral constraints will roughly be a function of the fluxes and colors of the covariant sources, along with the color prior itself. A scenario where color priors can help the model favor two sources is when merging two covariant sources results in a brighter source with colors less favored by the priors.

Using the source association procedure outlined in Section 4.1, one can calculate the mean error on position and flux for different two-source configurations. Figure 6 shows the mean absolute position errors for each source. The co-add and joint multiband runs recover more accurate source positions than corresponding single-band runs, which can be primarily attributed to higher S/N. Also plotted with dashed lines are the estimated position errors if one considers a source with no covariant neighbors and assumes that the likelihood dominates over the prior (see Appendix C of Portillo et al. 2017 for a derivation). As the separation between both sources increases, the position errors for single-band and multiband runs converge toward these lower bounds. This agreement demonstrates the robustness of probabilistic cataloging, especially when one considers that the sources are still highly covariant at these separations.

Zoom In Zoom Out Reset image size

Similar improvements can be seen in Figure 7, which shows the fractional flux error for both sources as a function of separation. Errors in flux may come from one sample source absorbing the flux from two highly covariant mock sources, or from two covariant sample sources with incorrectly distributed fluxes.

Zoom In Zoom Out Reset image size

SDSS is an imaging and spectroscopic survey designed to probe large-scale structure (York et al. 2000). Located at the Apache Point Observatory in New Mexico, the SDSS telescope is a 2.5 m, f/5 modified Ritchey-Chrétien wide-field altitude-azimuth telescope. Since it began collecting data in 2000, SDSS has covered over one-third of the sky, with over 500 million objects detected photometrically. Of those, nearly three million have observed spectroscopically with either the APOGEE spectrograph, used to study galactic structure in the Milky Way (Wilson et al. 2019), or the BOSS spectrograph, which gathers galaxy redshift information to probe cosmic expansion through baryon acoustic oscillations (Eisenstein et al. 2005).

5.1.1. SDSS Photometric Observation

5.1.2. SDSS Astrometric Calibration

In order to map pixel coordinates to celestial coordinates, the SDSS processing pipeline divides the total drift scanned data into frames and proceeds with astrometric calibration within each frame. For each frame, the pipeline first makes corrections to the original pixel coordinates (x, y) using a smooth cubic fit.

For color < (color)₀:

$$\begin{eqnarray}&&x^{\prime} =x+{g}_{0}+{g}_{1}y+{g}_{2}{y}^{2}+{g}_{3}{y}^{3}+{p}_{x}(\mathrm{color})\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&y^{\prime} =y+{h}_{0}+{h}_{1}y+{h}_{2}{y}^{2}+{h}_{3}{y}^{3}+{p}_{y}(\mathrm{color}).\end{eqnarray} \tag{ 19 }$$

For $\mathrm{color}\geqslant {(\mathrm{color})}_{0}$:

$$\begin{eqnarray}&&x^{\prime} =x+{g}_{0}+{g}_{1}y+{g}_{2}{y}^{2}+{g}_{3}{y}^{3}+{q}_{x}\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&y^{\prime} =y+{h}_{0}+{h}_{1}y+{h}_{2}{y}^{2}+{h}_{3}{y}^{3}+{q}_{y}.\end{eqnarray} \tag{ 21 }$$

The coefficients $\{{g}_{i},{h}_{i}{\}}_{i=1}^{3}$ and $\{{p}_{j},{q}_{j}{\}}_{j\in \{x,y\}}$ are fit for each frame.

These transformations from Pier et al. (2003) correct for optical distortions that are a function of column only. This is the direction of the drift scan and of differential chromatic refraction (DCR). The variation from optical distortions is typically of order ∼0.2 pixels, so the cubic fit accounts for that properly. DCR is modeled either as a linear function of color (r − i for r, i, and z; g − r for g; and u − g for u) or as a constant depending on the observing band and the color of the observed source. Once the pixel coordinates are corrected, pixels are mapped to Catalog Mean Place (CMP) celestial coordinates (μ, ν), the angular offsets along and perpendicular to the idealized great circle being scanned along. This is done through an affine transformation:

$$\begin{eqnarray}&&{\mu }_{\mathrm{CMP}}=a+{bx}^{\prime} +{cy}^{\prime} \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{\nu }_{\mathrm{CMP}}=d+{ex}^{\prime} +{fy}^{\prime} .\end{eqnarray} \tag{ 23 }$$

The mappings of these transformations are contained for each frame within SDSS asTrans files.⁸

Consider the pixel transformation between bands a and b. To transform ${{\mathbb{R}}}^{2}:{(x,y)}_{a}\to {(x,y)}_{b}$, we separate (x, y)_a into their integer and noninteger parts ${({x}_{0}+{dx},{y}_{0}+{dy})}_{a}$ by truncation and use the following:

$$\begin{eqnarray}&&{x}_{b}^{{\prime} }={x}_{b}({x}_{a,0},{y}_{a,0})+\left(\displaystyle \frac{{{dx}}_{b}}{{{dx}}_{a}}\right){{dx}}_{a}+\left(\displaystyle \frac{{{dx}}_{b}}{{{dy}}_{a}}\right){{dy}}_{a}\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&{y}_{b}^{{\prime} }={y}_{b}({x}_{a,0},{y}_{a,0})+\left(\displaystyle \frac{{{dy}}_{b}}{{{dx}}_{a}}\right){{dx}}_{a}+\left(\displaystyle \frac{{{dy}}_{b}}{{{dy}}_{a}}\right){{dy}}_{a}.\end{eqnarray} \tag{ 25 }$$

In the equations above, ${x}_{b}({x}_{a,0},{y}_{a,0})$ and ${y}_{b}({x}_{a,0},{y}_{a,0})$ are the integer-value transformed pixels taken directly from the asTrans files. The partial derivatives used in Equations (24) and (25) are approximated for each pixel with a first difference:⁹

$$\begin{eqnarray}&&\displaystyle \frac{{{dx}}_{b}}{{{dx}}_{a}}=\displaystyle \frac{{x}_{b}({x}_{a,0}+1,{y}_{a,0})-{x}_{b}({x}_{a,0}-1,{y}_{a,0})}{2}\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{dx}}_{b}}{{{dy}}_{a}}=\displaystyle \frac{{x}_{b}({x}_{a,0},{y}_{a,0}+1)-{x}_{b}({x}_{a,0},{y}_{a,0}-1)}{2},\end{eqnarray} \tag{ 27 }$$

where dx_b/dx_a and dx_b/dy_a are evaluated at $(x,y)=({x}_{a,0},{y}_{a,0})$.

This interpolation scheme allows us to map subpixel values by leveraging the information from local astrometric variation. In practice, these quantities are precomputed—six arrays specify ${[{x}_{0}]}_{b}$, ${[{y}_{0}]}_{b}$, $\left[\tfrac{{{dx}}_{b}}{{{dx}}_{a}}\right]$, $\left[\tfrac{{{dx}}_{b}}{{{dy}}_{a}}\right]$, $\left[\tfrac{{{dy}}_{b}}{{{dx}}_{a}}\right]$, and $\left[\tfrac{{{dy}}_{b}}{{{dy}}_{a}}\right]$ for pixels in the desired region. In this scheme a pixel transformation is computed by six array lookups followed by multiplication and addition operations, rather than the typical trigonometric transformations used in WCS evaluations.

Figure 8 shows the residual errors of linearized transformations from r to i and g, compared to the direct mappings from sky to i or g bands. The errors in x and y are centered about zero and below 10⁻⁴ pixels. As part of validation we also ran "round trip" tests, using separate linearized transformations to take a set of coordinates from $r\to g$ and then back from $g\to r$. The residuals in positions from that test were below 10⁻⁵ pixels and centered about zero, likely the result of numerical round-off errors. While there may be larger systematic errors that affect the astrometry, these tests validate our method as consistent and unbiased.

Zoom In Zoom Out Reset image size

Figure 9 shows a comparison of coordinate transformation times per source between astropy.wcs, a commonly used routine, and our linearized asTrans mappings. As the number of simultaneously perturbed sources increases, linearized transformations quickly outperform astropy.wcs; for N_src = 10³, there is a 40-fold speed improvement, and this grows to a factor of 80 by N_src = 10⁴. To emphasize the importance of having fast coordinate transformations, Figure 10 shows a breakdown of computational resources used in a run where ∼10³ sources were inferred. While the linearized coordinate transformations compose ∼5% of all computations, astropy.wcs would use ∼200% of the total computational resources used by PCAT in its current form.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this section, we present the results of multiband probabilistic cataloging applied to SDSS imaging of the globular cluster Messier 2 (M2). M2 is located 11.5 kpc away and contains about 150,000 stars with an overall stellar spectral type of F4 (Harris 1996, 2010 edition). Our tests focus on the same 100 × 100 pixel r-band cutout used in Portillo et al. (2017), with the addition of SDSS i- and g-band data. The cutout comes from run 2583, camcol 2, field 136, with pixel coordinate (x₀, y₀) = (310, 630) in the frame file defining the lower left corner. The center of the cutout in celestial coordinates is given by R.A. = ${21}^{{\rm{h}}}{33}^{{\rm{m}}}20\buildrel{\rm{s}}\over{.} 5$, decl. = $-0^\circ 48^{\prime} 26\buildrel{\prime\prime}\over{.} 5$ This field serves as a good test case because it is extremely crowded. It is so crowded, in fact, that the SDSS photometric pipeline Photo originally timed out while trying to catalog the region. An et al. (2008) produce a catalog of the same region using DAOPHOT on SDSS ugriz data. Sarajedini et al. (2007) provide a catalog of M2 derived from the HST Advanced Camera for Surveys (ACS). Both catalogs serve as important points of comparison. The DAOPHOT catalog of SDSS data represents results from a traditional approach to crowded field photometry. HST detects sources as faint as r ∼ 30 because of superior sensitivity and an angular resolution ∼20 times better than the Sloan telescope, meaning that its catalog can be treated as our reference "truth" catalog. Along with smaller pixel sizes, observations from HST are not affected by atmospheric seeing. The region we use has been observed by HST in the ACS F606W and F814W bands. We compare with F606W, which covers a similar bandpass to the SDSS r band.

The bias and gain of each band are obtained from the opECalib CCD electronics calibration file.¹⁰ Next, the asTrans coordinate mappings described in Section 5 are obtained for r → i and r → g. First, we generate mappings in which the mean color term is taken to be zero. We then incorporate the DCR correction by extracting coefficients p_x and p_y of Equation (18) from the asTrans files. The threshold colors ${(g-r)}_{0}\,=\,1.5$ and ${(r-i)}_{0}\gg 1$ are sufficiently large that we assume that all sample positions can be corrected with Equation (18). DCR corrections are computed for each source at every step in the Markov chain and are typically at the level of 10⁻² pixels.

In this work, a fixed PSF template is used for each band. This assumes that the PSF does not vary over the region of interest. Because the region we analyze is a 100 × 100 pixel region of M2 covering roughly $0.43\ {\mathrm{arcmin}}^{2}$, such an assumption is reasonable. The PSF templates in r, i, and g are extracted from SDSS psField files and are upsampled using a Lanczos kernel. As shown in Portillo et al. (2017), the quality of the PSF has a dramatic effect on the quality of the fit. This is important because the PSF fitting done by SDSS is unvalidated for crowded fields like the one we examine. While the SDSS pipeline does a good job fitting the r- and i-band PSFs, we find that the g-band PSF is notably broader, and inspection of residuals confirms that it is poorly estimated. To address this, we use the crowded field photometry code crowdsource (Schlafly et al. 2018) to refit the g-band PSF. We fit the PSF using a larger 600 × 600 pixel region centered on our image for more constraining power.

Catalog ensembles for two multiband runs are obtained, one using r- + i-band and one using r- + g-band data. To ensure that samples in the catalog ensemble are drawn from the posterior distribution, we discard the first 1500 of 3000 thinned samples of the chain as burn-in. In the case of r + i, PCAT infers 1380 ± 10 sources, while for r + g 1233 ± 5 sources are inferred. The corresponding single-band catalog ensemble contains ≈1100 sources in the same region. DAOPHOT only detects 356 sources in this region, whereas HST detects 6051 sources down to F606W ∼ 30, 1428 of which have F606W < 23.

In our first results, we observed oversplitting of some bright sources. We attribute this to astrometric miscalibration across bands—Pier et al. (2003) report relative astrometric accuracy of ∼25–35 mas between r band and u, g, i, z. As discussed in Section 3.2, astrometric miscalibration at the level of 10⁻² pixels is significant enough to affect sources at the bright end. To correct for this miscalibration, we take the brightest sources down to r = 19 and, for each source, absorb any neighboring sources within 1 pixel. Upon inspection, the oversplit adjacent sources are almost always low-significance sources. Once these are recombined, the corresponding fluxes and positions are recalculated. While we were able to reduce oversplitting in two-band runs with this correction, astrometric miscalibration had a more severe impact on the r + i + g three-band joint fit, the results of which we do not include in this work.

Color–magnitude diagrams for the two-band runs are shown in Figure 11. To obtain uncertainties for individual sources, we reduce the last 300 catalog samples to a "condensed catalog," which naturally marginalizes over the posterior catalog ensemble while providing a labeling that enumerates sources across different samples like a traditional catalog does. The process for producing a condensed catalog is outlined in Section 5 of Portillo et al. (2017). We institute a modest cut on the condensed catalog, removing sources present in less than 10% of samples. These are spurious sources that do not improve the catalog. Bright sources with large error bars are likely oversplit sources that are not converged. Also plotted are DAOPHOT catalog sources (green) and the fiducial sequence from An et al. (2008) for the full globular cluster.

Zoom In Zoom Out Reset image size

Posterior color distributions for each two-band run are displayed in Figure 12. While a Gaussian prior $\pi (r-i)\sim { \mathcal N }(0.25,1.0)$ is used, the posterior skews toward positive colors. Fainter sources in the catalog ensemble are less constrained by the data, which helps explains why brighter catalog sources have a posterior (black) that is narrower and more consistent with that of DAOPHOT catalog sources (green).

Zoom In Zoom Out Reset image size

To compare our catalog ensemble with the HST catalog, we associate sample sources with HST sources using two criteria: if a source is within 0.5 pixels of the HST position and has an r magnitude within 0.5 mag of the F606W magnitude, it is considered a match. These matching criteria are motivated by the fact that F606W is a wider band than r, and SDSS sources will be biased brighter owing to faint neighboring sources HST can resolve but SDSS cannot. In the context of source detection, the estimates of completeness and false discovery are not significantly affected by these criteria. If one assumes that the positions of sources in our region are generated from a spatial Poisson process (i.e., uniformly throughout the region), the expected fraction of spurious associations with the HST catalog given our criteria is ∼3% at r = 22 and drops to <1% for r < 20.

Figure 13 compares the completeness and false discovery rates of the DAOPHOT catalog (black) and single-band catalog ensemble from Portillo et al. (2017) (green). The single-band probabilistic catalog goes more than a magnitude deeper than DAOPHOT while maintaining a lower false discovery rate. Also plotted are the completeness/false discovery rates from the joint r + g (blue) and r + i (red) condensed catalogs. Multiband catalog fits further improve the completeness by several tenths of a magnitude while yielding lower false discovery rates. These results are roughly consistent with both the S/N estimates and multiband mock tests from Section 3. The improved deblending of the two-band condensed catalogs can be seen clearly between 18th and 20th magnitude. In this range, over 40% of all DAOPHOT catalog sources are false positives. The single-band probabilistic catalog from Portillo et al. (2017) achieves a false discovery rate (FDR) between 20% and 30% in this range, while our two-band condensed catalogs are consistently under 20%.

Zoom In Zoom Out Reset image size

Figure 14 shows the binned magnitude histograms for Hubble, PCAT, and DAOPHOT catalog sources. It is clear that on the faint end PCAT detects many more sources that match with Hubble catalog sources than DAOPHOT does. The same improvements in completeness and false discovery seen in Figure 13 are also reflected by these magnitude distributions.

Zoom In Zoom Out Reset image size

A useful quantity for comparing catalogs is the degradation factor, which is defined in Portillo et al. (2017) as the ratio of reported uncertainty in flux for a given source and expected uncertainty for an isolated source with the same flux. Sources with large degradation factors are those with large uncertainties, which are more likely to be heavily contaminated by neighbors or to be artifacts. As such, making successive cuts on the degradation factor is one way to describe the completeness–false discovery trade-off for our cataloger. Figure 15 shows completeness–false discovery rate curves for four magnitude bins between r = 17.5 and 21.5, with each curve reporting values for one- and two-band condensed catalogs. For the magnitude bin centered on r = 21, both r + i and r + g outperform the single-band condensed catalog—one can tell that these are better because their ROC curves are closer to the upper left portion of the completeness–false discovery plane, which is ideal. For r = 20, r + i is better than single-band, but the difference is smaller. Both multiband fits are worse at r = 19, and at r = 18 the r + g fit is worse. This decrease in performance may arise from the oversplitting of bright sources. Our estimates of completeness and false discovery come from one 100 × 100 pixel SDSS field of view, so the calculated ROC curves are also affected by sampling noise.

Zoom In Zoom Out Reset image size

This implementation of multiband probabilistic cataloging uses various computational and algorithmic optimizations to improve convergence. With half an hour of CPU time, it recovers similar results as the implementation from Portillo et al. (2017) does in 12 CPU-days, a speedup factor of over 500. This brings PCAT significantly closer to being a feasible pipeline for near-future surveys. Graphics processing units (GPUs) have been considered as a way to speed up probabilistic cataloging. However, many of the operations that would benefit from speed improvements are nontrivial to parallelize. For example, when sources are added or removed from an image, they need to be modified sequentially. There are modern CPUs that do this type of sequential operation efficiently. Another consideration is speed of matrix multiplication, which we use heavily in model evaluation (for subpixel shifts). Our matrix problem has a small internal dimension of 10, causing standard libraries (e.g., CUDA) to run at low GPU utilization. A modern Intel CPU (e.g., Sky Lake) with dual AVX-512 fused multiply-add engines can perform 64 floating-point operations per cycle per core. Our PCAT code nearly saturates that theoretical limit even with a small internal dimension. Because we have a much higher utilization on CPUs versus GPUs (10× or more), we currently obtain better performance on CPUs. Future versions of PCAT that allow for a spatially varying PSF might have a much larger internal dimension in the matrix multiplication and scale better on GPUs.

A crucial feature of any cataloger used for large-scale surveys is the ability to model galaxies, which are present in most fields in the sky. Galaxies with parameterized shapes (e.g., Sérsic profile) can be incorporated into PCAT's forward model; however, the higher dimensionality required to model galaxies makes sampling more challenging. Furthermore, such parameterizations may not be sufficient for galaxies with complicated morphologies. Future development of PCAT will focus on constructing efficient galaxy proposals.

The challenges of astronomical cataloging are tied to a number of more general problems within statistics. Through our prior that penalizes additional sources, we employ a form of model regularization, in which models are selected based on goodness of fit and model complexity. By sampling across models, we calculate relative model evidences (i.e., Bayes factors) by the relative abundance of samples from each model. Explicitly calculating the (normalized) model evidences is often computationally intractable in high-dimensional problems. Approximate methods exist such as variational inference (VI), which uses optimization techniques to maximize an evidence lower bound. VI excels at handling large data sets, where one wishes to quickly explore multiple models, but generally underestimates the variance of the posterior density (Blei et al. 2016; Regier et al. 2019).

Another alternative is to implement Hamiltonian reversible jump Monte Carlo, as has been done in Sen & Biswas (2017). In their framework, transdimensional moves have Hamiltonian Monte Carlo moves added at the beginning and/or end. The transdimensional move may then not require much tuning, as the Hamiltonian moves will guide the proposal to higher likelihoods. Preliminary tests of fixed-dimensional HMC on cataloging have been done with promising results.

Consider a catalog with N sources observed in k bands. The catalog may be represented as $C=\{(x,y,{f}_{1},{f}_{2},\ldots ,{f}_{k}){\}}_{n\,=\,1}^{N}$, as in Section 2.1. In the single-band case, one would simply include a prior π(f₁) on the flux of each source in that band. Trivially, one could extend this to multiple bands, imposing the same type of power law prior to each flux and using those in conjunction with one another. Instead of this, we impose a single reference band flux prior (r band in this work) along with color priors relating the remaining bands to the reference band. This takes the following form:

$$\begin{eqnarray}&&\pi ({\boldsymbol{f}})=\pi ({f}_{1})\times \prod _{i=2}^{k}\pi ({s}_{i}),\end{eqnarray} \tag{ 28 }$$

in which the color ${s}_{i}=-2.5{\mathrm{log}}_{10}({f}_{i}/{f}_{1})$. The color prior allows us to make use of the fact that the fluxes of a source observed in multiple bands are not independent. We choose a Gaussian prior on color $\pi ({s}_{i})\sim { \mathcal N }({\mu }_{i},{\sigma }_{i})$. The means μ_i and Gaussian widths σ_i depend on how well individual colors are constrained within certain astronomical populations. In practice, relatively wide Gaussian widths are used, only disfavoring sources with extremely atypical colors (e.g., $| r-i| \gt 4$). When new sources are created in birth moves, these color priors can be useful. Rather than draw independent flux samples for a new source, we choose to draw one flux f₁ and then choose other fluxes f_i by sampling from the colors π(s_i) and using those to calculate f_i. This results in more efficient proposals where birthed sources are drawn from a reasonable stellar locus.

To avoid overfitting our catalog model, PCAT institutes a prior on N, the number of sources. This parsimony prior, also referred to in the literature as a regularization prior, penalizes each additional source in the model. For an idealized Gaussian problem in one band, adding a source corresponds to an average improvement in log-likelihood of the maximum likelihood solution by 3/2 (1/2 per degree of freedom). This is a result of Wilks's theorem, which states that for a likelihood ratio Λ between hypotheses Θ and Θ₀, the quantity 2 logΛ will asymptotically be chi-squared distributed as χ² with degrees of freedom equal to the difference in dimensionality between Θ and Θ₀ (Wilks 1938). The parsimony prior that counteracts this overfitting effect is $\pi (N)\propto \exp \left(-\tfrac{3}{2}N\right)$. Because we add parameters to our model for each additional band, the prior is modified to

$$\begin{eqnarray}&&\pi (N)\propto \exp \left(-\displaystyle \frac{N}{2}(2+{n}_{\mathrm{bands}}\right).\end{eqnarray} \tag{ 29 }$$

While a Poisson prior on the mean number of sources N may be more consistent with our Bayesian hierarchical model, it is also less punishing of new sources. Tests using the Poisson prior produced a catalog ensemble with many sources that did not affect bright sources and had low S/N. If one specifically wanted to estimate the faint-end source population, the parsimony prior would not be appropriate. However, the focus of this work is deblending, and the parsimony prior in Equation (29) helps cut down on CPU time requirements because the penalty for new sources is stronger.

With these changes in effect, our prior is proportional to the following:

$$\begin{eqnarray}\begin{array}{rcl}\pi ({\boldsymbol{\theta }}) & \propto & \pi (N)\displaystyle \prod _{n=0}^{N}\pi ({x}_{n},{y}_{n})\pi ({{\boldsymbol{f}}}_{n})\\ & = & \pi (N)\displaystyle \prod _{n=0}^{N}\pi ({x}_{n},{y}_{n})\pi ({f}_{n,1})\displaystyle \prod _{i=2}^{k}\pi ({s}_{n,i}).\end{array}\end{eqnarray} \tag{ 30 }$$

In the case of multiple bands, we have some degree of freedom regarding how we can propose source fluxes. When proposing new sources, our implementation draws the reference band flux from a power law, while the remaining fluxes are calculated through fair draws on Gaussian color priors. To perturb the fluxes of existing sources, we draw proposals directly on all fluxes, such that our proposal distribution resembles the eventual stationary distribution on individual fluxes.

The size of position-shifting proposals also depends on source fluxes and is modified in the multiband case. While the optimal step size for position proposals is calculated using quantities relevant to the image, we make a modification that scales for more degrees of freedom per source. Fixing N_src somewhere on the order of N_max (2000 in this work),

$$\begin{eqnarray}&&{\sigma }_{x}\propto \displaystyle \frac{1}{{N}_{\mathrm{src}}(2+{n}_{\mathrm{bands}})}\displaystyle \frac{1}{{f}_{1}}.\end{eqnarray} \tag{ 31 }$$

To maintain detailed balance, ${f}_{1}=\max ({f}_{1,\mathrm{current}},{f}_{1,\mathrm{proposed}})$, which is conserved between a proposed step and its inverse.

In merge/split proposals, we would like to conserve the flux of our original source/sources. To some extent, we would like to conserve the color of our original source as well, in the sense that we would like to propose sources with reasonable colors. If the fraction of original source flux is drawn randomly as ${F}_{k}\sim \mathrm{Unif}[0,1]$ for each band separately, the two sources produced from a split are likely to have unreasonable colors, in which case our prior on color will suppress the proposal acceptance factor. Our proposal generates sources with colors perturbed from the original source color. In the reference band, we draw $\rho \sim \mathrm{Unif}[0,1]$ and calculate ${F}_{1}=({f}_{\min }/f)+\rho (1-2{f}_{\min }/f)$. This is designed so that as f approaches 2f_min (the lowest permissible to be a candidate for a split), F₁ goes to 0.5, and for $f\gg {f}_{\min }$, F₁ approaches ρ. For each band k other than the reference band, we generate a delta color ${\rm{\Delta }}s\sim { \mathcal N }(0,{\sigma }_{s})$, from which the flux fraction is calculated as

$$\begin{eqnarray}&&{F}_{k}=\displaystyle \frac{\exp \left(\tfrac{{\rm{\Delta }}s}{\kappa }\right){F}_{1}}{1-{F}_{1}+\exp \left(\tfrac{{\rm{\Delta }}s}{\kappa }\right){F}_{1}},\end{eqnarray} \tag{ 32 }$$

where $\kappa =\tfrac{2.5}{\mathrm{ln}10}\approx 1.08.$ One source receives fraction F_k of the original source flux in band k, while the other receives $1-{F}_{k}$. The colors of the two sources are then

$$\begin{eqnarray}&&{s}_{1,k}=\kappa \mathrm{ln}\displaystyle \frac{{F}_{k}}{{F}_{1}}\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&{s}_{2,k}=\kappa \mathrm{ln}\displaystyle \frac{1-{F}_{k}}{1-{F}_{1}}.\end{eqnarray} \tag{ 34 }$$

The acceptance factor for merges and splits also needs to be modified in the case of multiple bands. Consider the case of a proposal to split one source into two smaller sources. The acceptance ratio for a split can be expressed as

$$\begin{eqnarray}\begin{array}{rcl}{\alpha }_{\mathrm{split}} & = & \displaystyle \frac{{\pi }_{1}(x,y,{f}_{1},\{{s}_{i}{\}}_{i=2}^{n}){\pi }_{2}(x,y,{f}_{1},\{{s}_{i}{\}}_{i=2}^{n})}{{\pi }_{0}(x,y,{f}_{1},\{{s}_{i}{\}}_{i=2}^{n})}\\ & & \,\times \,\displaystyle \frac{{ \mathcal J }}{q({\rm{\Delta }}x,{\rm{\Delta }}y,\{{F}_{i}{\}}_{i=1}^{n})}\end{array}\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&=\ \displaystyle \frac{2\pi {k}^{2}}{A}\displaystyle \frac{{\pi }_{1}({f}_{1}){\pi }_{2}({f}_{1})}{{\pi }_{0}({f}_{1})q({F}_{1})}\prod _{i}\displaystyle \frac{{\pi }_{1}({s}_{i}){\pi }_{2}({s}_{i})}{{\pi }_{0}({s}_{i})q({\rm{\Delta }}{c}_{i})}{ \mathcal J }.\end{eqnarray} \tag{ 36 }$$

Here ${ \mathcal J }$ represents the Jacobian associated with our change of variables and q is the transition kernel from which proposals are drawn. The Jacobian in n bands is calculated to be

$$\begin{eqnarray}&&{ \mathcal J }=\prod _{i=1}^{n}\displaystyle \frac{\kappa }{{F}_{i}(1-{F}_{i})}.\end{eqnarray} \tag{ 37 }$$

The factor at the front of Equation (36) accounts for the kick range k at which the two split sources are separated, along with A, the area of the region being evaluated.

The corresponding acceptance factor associated with mergers is simply the reciprocal of the split likelihood factor, which maintains detailed balance across proposals. While the split/merge proposal used here is not fully optimized, it does maintain detailed balance and is sufficient for our purposes.

Floating background levels in a catalog fit can be difficult because the background is highly covariant with sources in a given image. In the case where the background is allowed to float as a free parameter in the catalog model, the optimal background proposal step size is given by ${\sigma }_{b}\sim {\sigma }_{\mathrm{pixel}}\sqrt{1/{N}_{\mathrm{pix}}}$ in units of photoelectrons. In this expression, ${\sigma }_{\mathrm{pixel}}$ is the pixel variance that depends on the data and gain, and N_pix is the total number of pixels in the region being perturbed. Most observational data sets give estimates of background for their images, though in the case of crowded fields these estimates are not always accurate because of the covariance mentioned above.

In cases when one wants to understand physical properties of an astrophysical background, floating the background level is an important element of the analysis. In the optical case, however, we are able to obtain a reasonable fit so long as the background is constrained at the ∼10% level. There is still covariance between background and astronomical sources in optical data, but exploring this covariance is not the focus of this analysis, and so for our purposes we choose to fix the background in our main runs. Future analysis of optical data with PCAT will likely consider this in more detail, especially because the background cannot be assumed to be spatially uniform over larger regions of sky.