Improved Point Source Detection in Crowded Fields using Probabilistic Cataloging

TL;DR

This paper introduces a probabilistic cataloging method for crowded astronomical fields, improving source detection accuracy and depth compared to traditional catalogs, especially in dense regions.

Contribution

It presents the first optical probabilistic catalog for crowded fields, capturing source covariance and deblending ambiguities, and introduces an algorithm to produce a condensed catalog similar to traditional ones.

Findings

Probabilistic cataloging recovers more sources than traditional methods.

It detects sources over a magnitude deeper than traditional catalogs.

The method achieves lower false discovery rates in crowded fields.

Abstract

Cataloging is challenging in crowded fields because sources are extremely covariant with their neighbors and blending makes even the number of sources ambiguous. We present the first optical probabilistic catalog, cataloging a crowded (~0.1 sources per pixel brighter than 22nd magnitude in F606W) Sloan Digital Sky Survey r band image from M2. Probabilistic cataloging returns an ensemble of catalogs inferred from the image and thus can capture source-source covariance and deblending ambiguities. By comparing to a traditional catalog of the same image and a Hubble Space Telescope catalog of the same region, we show that our catalog ensemble better recovers sources from the image. It goes more than a magnitude deeper than the traditional catalog while having a lower false discovery rate brighter than 20th magnitude. We also present an algorithm for reducing this catalog ensemble to a condensed catalog that is similar to a traditional catalog, except it explicitly marginalizes over source-source covariances and nuisance parameters. We show that this condensed catalog has a similar completeness and false discovery rate to the catalog ensemble. Future telescopes will be more sensitive, and thus more of their images will be crowded. Probabilistic cataloging performs better than existing software in crowded fields and so should be considered when creating photometric pipelines in the Large Synoptic Space Telescope era.

Figures

Equations

$$\theta =\{N,\{x_i,y_i,f_i{\}}_{i=1}^N,f_{max},I_{sky}\}.$$

$$\log \frac{\pi (N+1)}{\pi (N)}=-\frac{3}{2}.$$

$$\pi (N)\propto \exp \left(-\frac{3}{2}N\right).$$

$${\lambda }_{lm}=I_{sky}+\sum _{i=1}^N{f}_i\mathcal{P}(l-x_i,m-y_i).$$

$$\mathcal{L}=\prod _{l=1}^W\prod _{m=1}^H\frac{1}{\sqrt{2\pi {\lambda }_{lm}}}\exp \left(-\frac{(k_{lm}-{\lambda }_{lm}{)}^2}{2{\lambda }_{lm}}\right).$$

$$\log \mathcal{L}=\sum _{l=1}^W\sum _{m=1}^H-\frac{1}{2}\log {\lambda }_{lm}-\frac{(k_{lm}-{\lambda }_{lm}{)}^2}{2{\lambda }_{lm}}.$$

$$\log \mathcal{L}\approx \sum _{l=1}^W\sum _{m=1}^H-\frac{(k_{lm}-{\lambda }_{lm}{)}^2}{2{\lambda }_{lm}}.$$

$$FDR=\frac{FP}{FP+TP}$$

$${\frac{\partial \log \mathcal{L}}{\partial f}|}_{f=f_{\ast }}=0$$

$${\sigma }_f^2=⟨{\left(f-f_{\ast }\right)}^2⟩=\frac{{\int }_{-\infty }^{\infty }\exp \left(\log \mathcal{L}\right)(f-f_{\ast }{)}^{2}df}{{\int }_{-\infty }^{\infty }\exp \left(\log \mathcal{L}\right)df}$$

$${\sigma }_f^2\approx \frac{{\int }_{-\infty }^{\infty }\exp \left(\log \mathcal{L}{|}_{f=f_{\ast }}+\frac{1}{2}{\frac{{\partial }^2\log \mathcal{L}}{\partial f^2}|}_{f=f_{\ast }}(f-f_{\ast }{)}^2\right)(f-f_{\ast }{)}^{2}df}{{\int }_{-\infty }^{\infty }\exp \left(\log \mathcal{L}{|}_{f=f_{\ast }}+\frac{1}{2}{\frac{{\partial }^2\log \mathcal{L}}{\partial f^2}|}_{f=f_{\ast }}(f-f_{\ast }{)}^2\right)df}={\left(-{\frac{{\partial }^2\log \mathcal{L}}{\partial f^2}|}_{f=f_{\ast }}\right)}^{-1}$$

$$\frac{{\partial }^2\log \mathcal{L}}{\partial f^2}\approx -\sum _{l=1}^W\sum _{m=1}^H\frac{{\mathcal{P}}^2(l-x,m-y)}{{\lambda }_{lm}}{\left(1+\frac{k_{lm}-{\lambda }_{lm}}{{\lambda }_{lm}}\right)}^2$$

$$\frac{{\partial }^2\log \mathcal{L}}{\partial f^2}\approx -\sum _{l=1}^W\sum _{m=1}^H\frac{{\mathcal{P}}^2(l-x,m-y)}{f\mathcal{P}(l-x,m-y)+I_{sky}}$$

$${\sigma }_f^2\approx {\left(\sum _{l=1}^W\sum _{m=1}^H\frac{{\mathcal{P}}^2(l-x,m-y)}{f_{\ast }\mathcal{P}(l-x,m-y^{\ast })+I_{sky}}\right)}^{-1}$$

$$\frac{{\partial }^2\log \mathcal{L}}{\partial x^2}\approx -f_{\ast }^2\sum _{l=1}^W\sum _{m=1}^H\frac{1}{{\lambda }_{lm}}\left(\left[1+{\left(\frac{k_{lm}-{\lambda }_{lm}}{{\lambda }_{lm}}\right)}^2\right]{\left[\frac{\partial \mathcal{P}}{\partial x}(l-x,m-y_{\ast })\right]}^2-\frac{{\lambda }_{lm}}{f_{\ast }}\frac{k_{lm}-{\lambda }_{lm}}{{\lambda }_{lm}}\left[1+\frac{1}{2}\frac{k_{lm}-{\lambda }_{lm}}{{\lambda }_{lm}}\right]\frac{{\partial }^2\mathcal{P}}{\partial x^2}(l-x,m-y_{\ast })\right)$$

$${\sigma }_x^2\approx \frac{1}{f_{\ast }^2}{\left(\sum _{l=1}^W\sum _{m=1}^H\frac{1}{f_{\ast }\mathcal{P}(l-x_{\ast },m-y_{\ast })+I_{sky}}{\left[\frac{\partial \mathcal{P}}{\partial x}(l-x_{\ast },m-y_{\ast })\right]}^2\right)}^{-1}$$

Full text

catalogs — methods: statistical — methods: data analysis — globular clusters: individual (M2)

1 Introduction

While telescopes measure the incident light from the sky and generate images, observational astronomy is often done not with those images but with catalogs constructed from those images. In cataloging, emission sources are identified and their properties (e.g. flux, color, spatial extent) are measured. These catalogs can be used to study the population of sources itself: for example, the age and metallicity of a globular cluster can be measured from a color-magnitude diagram of its stars. The sources may also be used as markers or beacons, as in the Argonaut map (Green et al., 2015) which uses the flux and color of stars to map the 3D distribution of the intervening dust reddening and dimming the stars. Source catalogs are also used by scientists interested in only a certain class of source to identify candidates with specific properties for further followup.

When the sources are bright and well-separated, the catalog contains the information in the original image, and interpretation of uncertainty estimates is straightforward. However, in crowded fields, deblending becomes harder, and traditional catalogs cannot capture all of the information in the original image. For example, consider an emission peak that is best explained by two dim sources but is still consistent with a single bright source. A traditional catalog only allows a single enumeration of the sources in the image and thus must choose only one of these possibilities. This deblending ambiguity is not passed on to the downstream analyses that use the catalog, possibly causing these analyses to underestimate their uncertainties. Ambiguous sources may be flagged and cut by these analyses, but this cut potentially throws away information contained in the image. Furthermore, the inferred position and flux of a source will affect those inferred for its close neighbors. This covariance cannot be captured by a traditional catalog that simply lists the position, with uncertainty, for each source.

In addition, a consequence of the improved depth of future telescopes, such as the Large Synoptic Survey Telescope (LSST), is that most images will be in the crowded-field limit. Thus, current photometry algorithms will no longer suffice. For instance, in discussing potential photometry algorithms for the LSST pipeline, Becker et al. (2007) note that “no algorithms are able to meet the [Science Requirement Document] in [Point Spread Function] photometry – the numbers consistently fall short by a factor of 2-3.” Consequently, alternative photometry algorithms for cataloging crowded stellar fields and crowded fields in general are crucial to future surveys and to astrophysical research as a whole.

We propose probabilistic cataloging to improve the fidelity of catalogs in crowded fields. Instead of deriving a single catalog from an image, probabilistic cataloging infers the posterior distribution of catalogs in a Bayesian framework. This posterior distribution is sampled repeatedly (each sample being a catalog) to yield a catalog ensemble. Each catalog in the ensemble is a fair draw from the posterior distribution and is thus consistent with the data, and the catalog ensemble represents the state of knowledge about what sources could be in the image. In the case of ambiguous deblending detailed previously, most catalogs in the ensemble would explain the peak of emission with two dim sources but some catalogs would contain just a single bright source, reflecting the fact that a single source is still an allowed possibility. The relative prevalence of the two-source and one-source explanations in the ensemble would reflect the statistical favorability of these explanations relative to each other. The derived catalog ensemble reflects this ambiguity, allowing propagation of deblending uncertainties to downstream analyses. Source-source covariance is naturally captured by a catalog ensemble, and resulting uncertainties are naturally marginalized over the effects of neighbors, including neighbors too faint to appear in a traditional catalog. This marginalization over faint sources is a non-trivial benefit. While not all users of the catalog will need a list of possible 4-sigma sources, many users will be interested in 20-sigma sources that have errors correctly marginalized over the effects of neighboring 4-sigma sources.

Several recent efforts have begun to implement probabilistic cataloging in astronomy. Hobson & McLachlan (2003) present it as “simultaneous detection of all objects” and apply it to a toy problem: a $200\times 200$ pixel image with 8 discrete 2D Gaussian shaped objects of varying flux and radius. Brewer et al. (2013) extend its application to a crowded field, applying it to a simulated $100\times 100$ optical image with 1,000 stars and showing that it outperforms SExtractor (Bertin & Arnouts, 1996). Jones et al. (2015) implement probabilistic cataloging on a $25"\times 25"$ Chandra field with 14 sources, incorporating spectral models for the sources in order to use photon energy information to help deblend sources.

Other Bayesian source detection algorithms have been implemented; unlike probabilistic cataloging, they do not produce a catalog ensemble. For example, Hobson & McLachlan (2003) present an iterative object detection method which stops when the Bayesian evidence disfavors adding additional sources. Carvalho et al. (2009) speed up this approach by a factor of $\sim 100$, and further refine it in Carvalho et al. (2012) to produce the Planck Early Release Compact Source Catalogue. These implementations are considerably faster than probabilistic cataloging and work well on uncrowded fields. However, these methods do not capture source-source covariances and make approximations that do not hold in crowded fields. Masias et al. (2012) offer a review of source detection approaches, including Bayesian approaches.

In this work, we present a probabilistic catalog inferred from part of a Sloan Digital Sky Survey (SDSS) (Gunn et al., 2006; Eisenstein et al., 2011; Ahn et al., 2012) image of the globular cluster M2, the first probabilistic catalog constructed from optical data. We evaluate both our probabilistic catalog and a traditional catalog derived from the same SDSS image by comparing to a much deeper catalog from the Hubble Space Telescope (HST). We also introduce a procedure to distill the probabilistic catalog samples into a condensed catalog which is easier to use in downstream analyses. Like a traditional catalog, it lists sources with their properties and uncertainties, but unlike a traditional catalog, the condensed catalog also marginalizes the sources’ properties over their neighbors and includes a measure of confidence in each source’s existence.

The rest of this paper is organized as follows. In Section 2, we describe the SDSS observations used to construct our probabilistic catalog, the traditional catalog we compare against, and the HST catalog we use as ground truth to assess both catalogs. In Section 3, we detail our Bayesian point source model and the MCMC sampler that we use. In Section 4, we present our probabilistic catalog and compare it against a traditional catalog. In Section 5, we present our condensed catalog and compare it to the full probabilistic catalog. We also compare the reported flux errors to the flux errors expected in a sparse field in order to identify spurious sources. In Section 6, we discuss implications and conclude.

2 Observations

In order to test the performance of probabilistic cataloging against traditional cataloging, we chose an image of the globular cluster M2 from SDSS. Globular clusters are rich with stars, providing challenging crowded fields to catalog. In addition, in globular clusters, nearly every source is a star rather than a galaxy, allowing the image to be modeled using only point sources. An et al. (2008) cataloged the same SDSS image of M2 using DAOPHOT (Stetson, 1987, 1994), providing a traditional catalog (referred to hereafter as the “DAOPHOT catalog”) to which we can directly compare our probabilistic catalog. In addition, the ACS Globular Cluster survey (Sarajedini et al., 2007) cataloged M2 using imaging from the Advanced Camera for Surveys on HST. HST has $\approx 20$ times the angular resolution of the Sloan Foundation 2.5 m Telescope, so the stars are much more widely separated in the HST image and are thus easily picked up by traditional cataloging methods. In addition, the HST image has $\approx 30$ times the exposure of the SDSS image (five 340 s HST exposures vs the 34.9 s SDSS exposure). We take the catalog made from the HST data (referred to hereafter as the “HST catalog”) as a reliable reference against which catalogs made from the SDSS data can be judged.

M2 is a rich and compact globular cluster containing about 150,000 stars located 11.5 kpc away. It has a core radius of 0.34 arcminutes (1.1 parsecs) and a half-light radius of 1.08 arcminutes (3.6 parsecs) (Harris, 1996, 2010 edition).

Our main test image is SDSS run 2583, field 136, camcol 2, in the r band. There is no SDSS catalog for this field, because the survey photometric pipeline, Photo (Lupton et al., 2001), timed out. For most of our tests, we focus on a $100\times 100$ pixel cutout from the image $\approx 2$ arcminutes from the center of the cluster. This patch is crowded, but away from the heavily saturated core. The HST catalog identifies 1,049 sources brighter than 22nd magnitude in the F606W band in this patch. F606W is much broader than SDSS r, but centered at roughly the same wavelength, and F606W magnitudes for main sequence stars agree with SDSS r magnitudes to a few tenths of a magnitude.

3 Methods

3.1 Model Space and Priors

The model space is the space of possible point-source catalogs. We assume that there are an unknown number, $N$, of point sources in the image. In our case of point sources in a single band, each source $i$ is described by its position $(x_{i},y_{i})$ and flux $f_{i}$. We assume the sources are uniformly distributed spatially and that their fluxes follow a power-law distribution between $f_{min}$, which is fixed, and $f_{max}$, which is unknown. The reciprocal of $f_{max}$ is given a uniform prior on $(0,f_{bound}^{-1})$, where $f_{bound}$ is slightly dimmer than the brightest point source in the image, forcing $f_{max}$ to be at least as bright as the brightest point source.

In addition to the point source emission, the image has unknown uniform sky emission $I_{sky}$, which is given a log-uniform prior. Thus, a single catalog can be described by a parameter vector $\theta$:

[TABLE]

Note that the catalog, as defined above, includes not only the point sources and their properties (which would constitute a “catalog” in normal usage), but also the hyperparameters describing the point source population and the sky level as a nuisance parameter. Also, a traditional catalog would include uncertainties on the positions and fluxes, but each sample of a probabilistic catalog does not have uncertainties on its parameters. Instead, the uncertainty manifests in the distribution of these parameters across the catalog ensemble.

The space of possible catalogs is transdimensional, that is to say, not of fixed dimension. With 3 parameters for each source and 2 other parameters, the space of possible catalogs with $N$ sources is $3N+2$ dimensional. The space of all possible catalogs can thus be thought of as the union of the spaces of possible catalogs with $N=0,1,2,3$, and so on. It thus contains fixed-dimensional subspaces each with a different dimensionality. Meaningful priors can still be assigned on this space while the likelihood is in fixed-dimension data space, so a posterior over the space of possible catalogs is definable. Such a transdimensional space can be sampled with methods like Reversible Jump MCMC (Green, 1995) or Birth Death MCMC (Stephens, 2000).

We treat the number of point sources, $N$, as arising from a Poisson process which would yield $\mu$ point sources on average in the imaged region. Daylan et al. (2017) is concerned with constraining populations of sub-threshold sources and treats $N$ as the answer to the question: “What is the number of point sources above a given minimum flux, $f_{min}$ (counting even sub-threshold sources)?” Thus, they place a log uniform prior on $\mu$, not wanting to penalize possible sub-threshold sources. However, the current work focuses on deblending overlapping and significant sources. We are not interested in the population of sub-threshold sources, and they will not greatly affect the deblending of highly significant sources. Since the time taken to evaluate the likelihood for each sample increases with the number of included point sources, we instead add a prior on $N$ that penalizes additional point sources and cuts out these sub-threshold sources. For a Gaussian problem, when all parameters are away from the boundaries and the data are not overfit, adding a source improves the log likelihood of the maximum likelihood solution, on average, by $\frac{3}{2}$ ($\frac{1}{2}$ for each of three parameters $x_{i},y_{i},f_{i}$). Thus we choose a prior that counteracts the expected log likelihood gain from additional sources:

[TABLE]

This choice implies an exponential prior on $N$:

[TABLE]

We set $f_{min}$ to correspond to the SDSS 95% completeness limit in a sparse field ($\approx 22$ mag in r band), expecting the dimmest significant source in the crowded field to be brighter than this limit.

3.2 Generative Model and Likelihood

To calculate the likelihood of a given catalog, the corresponding model image must be produced. The data are photoelectron counts $k_{lm}$ in a $W\times H$ grid of pixels with coordinates $(x,y)=(l,m)$. We use the pixel-convolved point spread function extracted by the standard SDSS pipeline for the center of our image, $\mathcal{P}(\Delta x,\Delta y)$, to predict the expected counts $\lambda_{lm}$ for each pixel from a uniform sky background $I_{sky}$ and nearby sources:

[TABLE]

In the limit of many photoelectron counts, the noise when $\lambda_{lm}$ photoelectrons are expected is Gaussian with standard deviation $\sqrt{\lambda_{lm}}$ photoelectrons. The likelihood is then

[TABLE]

Taking the log likelihood and dropping additive constants because MCMC sampling only requires likelihood ratios, $\log\mathcal{L}$ gives:

[TABLE]

The first term, being a logarithm, is much less sensitive to changes in the $\lambda_{ij}$ than the second term. Assuming the first term is constant reduces the log likelihood evaluations to only additions and multiplications, greatly speeding up sampling:

[TABLE]

3.3 MCMC Sampling

This work performs MCMC sampling of the posterior distribution using Diffusive Nested Sampling111We use the DNest3 C++ implementation, available at https://github.com/eggplantbren/DNest3. (Brewer et al., 2011). Rather than directly sampling the posterior distribution as in Metropolis-Hastings, Diffusive Nested Sampling constructs an alternative target distribution. The Diffusive Nested Sampling algorithm first samples the prior and then uses these samples to determine a likelihood threshold that contains $\approx e^{-1}$ of the prior mass above it. The sampler then samples from a modified version of the prior where the parameter space above the likelihood threshold (the first level) is upweighted. Using these samples, a second likelihood threshold is determined such that the parameter space above it (the second level) contains $\approx e^{-2}$ of the prior mass above it. The sampler then samples from a modified version of the prior where the first level is still upweighted, but the second level is upweighted more. This process continues until the next level constructed is determined to hold so little of the posterior mass that it is unimportant.

In constructing these levels, the sampler behaves similarly to simulated annealing, but with a cooling schedule tailored to the posterior distribution using samples. Once enough levels have been constructed, the sampler continues to sample from this series of levels. The resulting samples can later be reweighted to yield samples from the posterior. This target distribution allows the sampler to escape posterior peaks by sampling from a lower (and thus less constraining) level; the sampler is even allowed to return to the prior. This freedom allows the sampler to avoid being trapped in a single posterior peak and instead find multiple peaks in a multi-modal distribution.

4 Probabilistic Catalog

4.1 Computational Requirements

All code was run on a machine with two six-core Intel Xeon E5-2667 processors at 2.9 GHz. Running our probabilistic cataloging code using all 12 cores on the selected $100\times 100$ pixel image, inferring $\approx 1,100$ sources, required a day of wall clock time to yield useful results. Each sample in the MCMC chain requires 1 core-ms on average.

4.2 Representation as Catalog Ensemble

To clarify the representation of the probabilistic catalog as a catalog ensemble, we present fair samples from the ensemble in Figure 1. Sources that are inferred confidently appear in all of the sample catalogs, while less confident inferences only appear in a fraction of the samples. The uncertainty in the position and flux of a confidently identified source manifests as the distribution of these properties across sample catalogs. The uncertainty in deblending ambiguous sources is reflected in the differing deblendings made in sample catalogs.

4.3 Completeness

We take the HST catalog as having the true identities and positions of the stars in the selected patch of M2. Because of HST’s superior angular resolution, the HST image is much less crowded (in terms of sources per resolution element) than the SDSS image, so identifying stars is straightforward, and their positions are measured well. Although the formal error on the HST positions is negligible, we are still limited by how well the HST and SDSS astrometric solutions are matched, and there may be noticeable movement of high proper-motion stars.

In Figure 2, we present a comparison of the completeness of our catalog ensemble and the DAOPHOT catalog against the stars in the HST catalog. Completeness is plotted against HST 606W magnitudes, the closest HST band available in this field to SDSS r band. A source is taken to be a match if it is within 0.75 pixels of the HST position (cp. the full width at half maximum of the SDSS point spread function, 2.2 pixels) and its r magnitude is within 0.5 mag of the 606W magnitude. This matching criterion is loose for two reasons: first, 606W is a wider band than r; and second, some stars have faint neighbors that are detectable in the HST image but are completely blended in the SDSS image. Because the different samples of the catalog ensemble can contain differing sources, the catalog ensemble completeness is reported as the average of the completeness of the sample catalogs. If the ensemble has $M$ samples, this averaging is equivalent to counting a match in a single sample as $M^{-1}$ of a match. The DAOPHOT catalog starts losing completeness at 18th magnitude, becoming incomplete by 22nd magnitude. The catalog ensemble, however, goes over 1 magnitude deeper at every level of completeness.

4.4 False Discovery Rate

A catalog with higher completeness is not necessarily a better catalog unless it also has a comparable or lower false discovery rate. Phrased differently, a catalog with sources at all possible positions and fluxes trivially has 100% completeness, but is a useless catalog because nearly every source is a false positive. Borrowing terminology from the binary classification literature, we define the false discovery rate ($FDR$) in terms of the number of true positives ($TP$) and false positives ($FP$):

[TABLE]

with a true positive being a catalog source with an HST source satisfying the same match criteria as for completeness in Section 4.3, and a false positive being a catalog source without any such match. In other words, the FDR is the fraction of the sources in the catalog which are false positives. Similarly to when we calculate completeness, when considering $M$ samples from the catalog ensemble, we count each source in a sample catalog as $M^{-1}$ of a (true or false) positive.

In Figure 3, we present a comparison of the false discovery rate of the catalog ensemble and the DAOPHOT catalog, as a function of SDSS r band magnitude. The DAOPHOT catalog false discovery rate peaks at 47% at 18th magnitude and stays above 10% for fainter magnitudes. The catalog ensemble false discovery rate is lower than the DAOPHOT false discovery rate until 20th magnitude. The false positive rate rises because closely-separated pairs of dim sources are fit with a single brighter source which does not meet the match criterion for either dim source. Beyond 21nd magnitude, the DAOPHOT catalog cannot be assigned a false discovery rate because it does not identify any sources that faint.

4.5 Residuals

It is instructive to look more closely at some very crowded regions of the image. In Figures 4, 5, 6, and 7, we overplot 300 stacked probabilistic catalog samples on a $10\times 10$ pixel region of the SDSS image, along with the DAOPHOT catalog and the HST catalog. Consider Figure 4: the DAOPHOT catalog claims that there are five sources. Comparing to the HST catalog, three of the DAOPHOT sources correspond to bright HST sources, while two of them correspond to the blending of several closely spaced HST sources. The probabilistic catalog captures the same three bright HST sources, better deblends the closely spaced stars into five, and places additional sources corresponding to single HST sources or collections of closely spaced HST sources. The probabilistic catalog residuals are smaller than those of the DAOPHOT catalog. Because of the probabilistic catalog’s better deblending, the probabilistic catalog residuals are much smaller around crowded stars than the DAOPHOT catalog residuals.

5 Recovering a Condensed Catalog

The burden of retaining and processing a large number of catalog samples may be a serious obstacle to using probabilistic cataloging. There is little advantage to a catalog ensemble if sources are unambiguously detected, but even in cases with severe crowding, a large fraction of the sources may be unambiguous. In this case, it may be useful to reduce the catalog ensemble to a single catalog, which we term a condensed catalog. This catalog has the convenience of a traditional catalog of one entry per object with uncertainties on each parameter, but retains the advantage of a fully marginalized posterior. That is, the flux and position uncertainties account for the effects of close neighbors. This middle road between traditional methods and a full catalog ensemble may be advantageous for many applications, especially if the user is most interested in properties of the brightest sources properly marginalized over the faint sources. In this section we propose a way to recover a condensed catalog from a catalog ensemble, and assess the condensed catalog’s performance.

5.1 Labeling Sources

The main complication in reducing a catalog ensemble to a condensed catalog is that the probabilistic catalog’s prior and likelihood are invariant under relabeling of the sources. While the MCMC sampling works by perturbing the properties of an ordered list of sources, the ordering of sources in this list cannot be taken as a consistent labeling of sources across samples in the catalog ensemble. It is useful to assign labels across samples to sources that are confidently detected in the catalog ensemble (i.e. they appear with well-constrained properties in a significant fraction of the samples). This labeling allows the position and flux distribution of these sources to be reported, just as in a traditional catalog. With a labeling procedure that at least works for unambiguous sources, we can evaluate the catalog ensemble’s false discovery rate for the sources it infers confidently. Labeling becomes more difficult for faint and crowded sources; however, the catalog ensemble makes less confident inferences about these sources, and its false discovery rate for low-confidence inferences is less problematic. We present a simple labeling procedure that recovers the sources that the probabilistic catalog confidently infers.

Our labeling procedure starts by constructing a seed catalog that denotes likely positions for identified by the catalog ensemble. First, the sources from all of the samples in the catalog ensemble are stacked. A given source is considered a neighbor to the closest source from each of the other samples in the ensemble within 0.75 pixels, if such a source exists for that other sample. Then, the source with the most neighbors is denoted a seed (if there is a tie, a source is chosen arbitrarily). The new seed and its neighbors are removed from the seeding algorithm, preventing them from becoming future seeds and decrementing the neighbor count of their neighbors. The procedure of identifying seeds and removing them and their neighbors is repeated until all sources have been exhausted.

Then, we consider the brightest source in the seed catalog and identify the brightest matched source from each of the probabilistic catalog samples within a search radius (0.75 pixels), if there is such a source. Each match is labeled as corresponding to the brightest seed catalog source. These matched sources are eliminated from the matching algorithm, as they now have labels. This matching is then repeated for the second brightest source in the seed catalog, then the third brightest, etc. until all sources in the seed catalog have been considered. Using these labels, a condensed catalog can be produced containing a mean position and flux for each source, with error bars, as well as a prevalence (the fraction of samples containing this source).

5.2 Comparing to Catalog Ensemble

The condensed catalog naturally marginalizes each labeled source’s properties over the properties of nearby sources. However, in reporting independent distributions for each source, the condensed catalog contains less information than does the catalog ensemble. The condensed catalog does not capture source-source covariance like the catalog ensemble does, but it does reflect the impact that source-source covariance has on the uncertainty in each source’s properties. We investigate the impact of this information loss on the completeness of the condensed catalog against the HST catalog. We generalize the notion of completeness to condensed catalogs. We do so by counting a condensed catalog source with a prevalence $p$ that matches with an HST source not as one match but as a fraction $p$ of a match. If, for all HST sources, all of the catalog ensemble sources that matched that HST source carry the same label and no other sources carry this label, then this completeness will match the catalog ensemble completeness. In Figure 9, we present a comparison of the completeness of the probabilistic catalog, condensed catalog, and the DAOPHOT catalog. The condensed catalog performs almost identically to the probabilistic catalog, demonstrating that the labeling is stable for at least the sources that correspond to HST sources.

We also characterize the condensed catalog’s false discovery rate as compared to the catalog ensemble. To generalize the false discovery rate (Equation 8) to condensed catalogs, we count a condensed catalog source with a prevalence $p$ as a fraction $p$ of a true positive if it matches an HST source, or as $p$ of a false positive if it does not. If, for all HST sources, all of the catalog ensemble sources that matched that HST source carry the same label and no other sources carry this label, then this false discovery rate will match the catalog ensemble false discovery rate. The false discovery rate of the condensed catalog is compared to that of the catalog ensemble and DAOPHOT catalog in Figure 10. The condensed catalog’s false discovery rate is comparable to that of the catalog ensemble’s, showing that labeling does not greatly affect the false discovery rate.

5.3 Completeness—False Discovery Tradeoff

To help identify spurious sources and sources heavily contaminated by their neighbors, we compare the errors reported in the condensed catalog to the errors that would be expected for a given source based on its flux in a sparse field (see Appendix C for the derivation). Sources that are contaminated by their neighbors will have larger errors than they would in a sparse field. Additionally, sources that are not well fit by our model because they are unresolved blends or artifacts of PSF mismatch will also have larger errors than expected. Cutting sources with a high ratio of reported to expected error (denoted here as the “degradation factor” and abbreviated as DF) may then remove many spurious sources while losing few legitimate sources. With a very loose cut at a flux degradation factor of 8, the false positive rate falls dramatically between 17th and 20th magnitude while only modestly decreasing completeness. The flux degradation factor can also be calculated for the DAOPHOT catalog based on the reported magnitude errors. Note that the DAOPHOT flux errors are not marginalized over neighboring sources, so DAOPHOT sources generally have lower degradation factors and are thus less affected by degradation factor cuts. Making the same cut at degradation factor 8 decreases the false positive rate of the DAOPHOT catalog without visibly impacting the completeness. Still, the DAOPHOT catalog’s false discovery rate is higher than that of the condensed catalog with the same degradation factor cut. We plot the completeness for the condensed catalog and DAOPHOT catalog, with and without this degradation factor cut, in Figure 11, and we plot the false discovery rates in Figure 12.

Cuts on parameters like the degradation factor can be thought of as tuning parameters for a catalog. More permissive cuts allow more real sources to be identified, but at the risk of allowing more false positives. Conversely, more restrictive cuts decrease the number of false positives, but at the cost of losing completeness. Inspired by receiver operating characteristic curves in the binary classification literature, we plot the performance of the DAOPHOT and condensed catalogs in the space of completeness versus false discovery rate. Varying the degradation factor cut for a given catalog, and thus making a different trade off between high completeness and low false discovery rate, traces a curve in this space. A perfect cataloger would lie at the top left —full completeness with no false positives. A cataloger whose curve lies completely above and to the left of another cataloger’s curve is superior in both metrics. We report the completeness versus false discovery rate curves at 18th, 19th, 20th, and 21st magnitudes in Figure 13. The condensed catalog clearly outperforms the DAOPHOT catalog, except at 21st magnitude if a false discovery rate of less than 15% is demanded. However, the false discovery rate is not well-measured for these restrictive cuts that do not admit many source candidates: the DAOPHOT catalog false discovery rate is 6/40 compared to the condensed catalog’s 5/35. Other parameters could be used in addition to attempt to further improve the condensed catalog’s performance. We defer a full study of the completeness-false discovery tradeoff obtained with parameters other than flux degradation factor for future work.

6 Discussion and Conclusion

We implement a probabilistic cataloger that produces an ensemble of point source catalogs of a crowded SDSS field in M2. We demonstrate that this catalog ensemble better recovers sources than the DAOPHOT catalog does by comparing to the sources identified by the HST catalog in the same field. The catalog ensemble is more than 1 magnitude deeper in completeness and has a lower false discovery rate above 20th magnitude. Many closely spaced sources that are blended in the DAOPHOT catalog are deblended in the catalog ensemble, and the residuals of the catalog ensemble model image are smaller and less spatially structured than those of the DAOPHOT model image. DAOPHOT seeds candidate source locations using a peak finder on a PSF-convolved image, making it likely to blend closely spaced sources. By instead proposing new sources’ locations and evaluating their likelihood, probabilistic cataloging is better able to separate closely spaced sources. In a crowded field, better deblending is a great advantage, explaining the catalog ensemble’s superior performance in completeness, false discovery rate, and residuals.

Distilling the obtained catalog ensemble to a condensed catalog does not lose much of the information contained in the catalog ensemble about confidently inferred sources. The condensed catalog has comparable completeness and false discovery rate to the catalog ensemble. We define a metric of quality: the degradation factor, the ratio of the flux error to the expected flux error for a source if it were isolated. A loose degradation factor cut drastically reduces the false positive rate while not impacting completeness by removing sources with abnormally high error, including some artifacts of PSF mismatch and incorrectly blended sources. While the condensed catalog appears similar to a traditional catalog, the errors for a source in this condensed catalog have been marginalized over the possible properties of its neighbors and nuisance parameters like the sky level.

The labeling algorithm for the condensed catalog detailed in Section 5.1 works well because most of the sources inferred were confidently inferred. In choosing our prior on the number of sources, we filtered out less significant sources that would have been included with less confidence. Labeling probabilistic catalogs including these less confident sources will be more difficult.

While our current implementation of probabilistic cataloging is only applicable to point sources, it can be generalized to extended sources like galaxies. In addition to position and flux, each source would have shape parameters (e.g. the parameters of an exponential or de Vaucouleurs profile) which would be used to calculate the model image against which the data are compared. Star-galaxy separation for galaxies that are not much bigger than the PSF may require special attention.

The main disadvantage of probabilistic cataloging is that it takes much more computing time than traditional cataloging. Given the amount of time needed to process the chosen $100\times 100$ image, cataloging an optical sky survey would be impractical with the current implementation of probabilistic cataloging. However, the current implementation has not been optimized for computational speed, so substantial gains in speed can be expected as the implementation is further developed. Also, the proposals used in the MCMC are unoptimized and uninformed - studying how to make better proposals may significantly reduce the number of steps needed in the MCMC chain. Alternatively, there may exist some approximation of probabilistic cataloging that is much faster while still retaining its advantages in crowded fields.

As telescopes become increasingly sensitive and detect more sources, more of the sky will become crowded enough to pose problems for traditional cataloging algorithms. If probabilistic cataloging or some derivative of it can be sped up, near-future computing power may make it worth running on the most crowded fields of an optical survey, or perhaps the entire survey itself. In demonstrating that probabilistic cataloging has superior performance in crowded fields, we contend that it should be further developed in order to create a photometric pipeline that will be able to meet the challenges posed by the next generation of optical surveys.

SKNP is supported in part by a Natural Sciences and Engineering Research Council of Canada Postgraduate Scholarship. BCGL was supported in part by the Harvard College Research Program. We thank Brendon Brewer, Charles-Antoine Collins Fekete, Dan Foreman-Mackey, Gregory Green, Kiefer Hicks, Albert Lee, and Zachary Slepian for useful discussions during the course of the project. The computations in this paper were run on the Odyssey cluster supported by the FAS Division of Science, Research Computing Group at Harvard University. This research has made use of NASA’s Astrophysics Data System. Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.S. Department of Energy Office of Science. The SDSS-III web site is http://www.sdss3.org/. SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, University of Cambridge, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofisica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University.

Appendix A MCMC Proposals

Choose a number of sources to perturb, and then perturb their positions and fluxes. 2. 2.

Perturb flux distribution parameters, changing all source fluxes to remain at the same quantile. 3. 3.

Perturb flux distribution parameters, keeping all source fluxes fixed. 4. 4.

Choose a number of sources to add, and then add sources drawn from the prior with the current flux distribution parameters. 5. 5.

Choose a number of sources to remove, and then remove sources chosen at random.

Whenever a number of sources is chosen, a multi-scale distribution is used with a minimum of 1 and a maximum of order the maximum number of sources allowed (3,000 in this work). Whenever a parameter is perturbed, a multi-scale distribution is used with a smallest step of order $10^{-6}$ the prior width and a largest step of order the prior width.

Appendix B Sparse Field Validation

To validate our photometry, we also run our probabilistic cataloger on a sparse field where SDSS Photo performs well. A significant fraction of sources in this sparse field are galaxies, and thus extended. Our cataloger assumes all sources are point sources, so we only compare photometry for sources that Photo identifies as point sources. We find all of the point sources that Photo finds, and our photometry agrees with Photo to 0.015 mag at the bright end, growing to 0.047 mag by 20th magnitude (see Figure 14). The positions that we infer agree with the Photo positions to 0.012 pixels at the bright end and 0.058 pixels by 20th magnitude (see Figure 15).

Appendix C Degradation Factor

To simplify the derivation of the expected flux error, we will assume that the likelihood dominates over the prior and that the approximations made in Equation 7 hold. Assuming that the position of this source, the positions and fluxes of the other sources, and the nuisance parameters are well known, the maximum likelihood flux for a single source $f_{*}$ can be found by maximizing the log likelihood:

[TABLE]

The flux variance about this maximum likelihood value is then:

[TABLE]

Taylor expanding the log likelihood about $f_{*}$ to second order and using Equation C1 yields:

[TABLE]

Assuming this source is well isolated, the model flux is $\lambda_{lm}\approx f\mathcal{P}(l-x,m-y)+I_{sky}$. Taking the second derivative of Equation 7 yields:

[TABLE]

Assuming the fractional deviations $(k_{lm}-\lambda_{lm})/\lambda_{lm}$ to be small gives:

[TABLE]

Substituting this expression into Equation C3 gives

[TABLE]

The reported errors in the condensed catalog are plotted in Figure 16. Adding a noise floor of 1% to account for effects like PSF mismatch which this derivation ignores, the expected flux error gives an excellent lower bound to the reported flux errors. Some sources near $f_{min}$ appear to have smaller flux errors because their flux distributions are cut off at the faint end by $f_{min}$. Because of how crowded the field is, many sources have flux errors many times larger than expected in the sparse limit.

By a similar procedure, the position uncertainty can be found by first taking the second derivative of Equation 7 with respect to $x$ or $y$:

[TABLE]

The second term is suppressed by the pixel model flux $\lambda_{lm}$ divided by the entire flux of the star $f_{*}$. Dropping this term and assuming the fractional deviations $(k_{lm}-\lambda_{lm})/\lambda_{lm}$ to be small gives:

[TABLE]

Adding a noise floor of 0.005 pixels in both $x$ and $y$, the expected position error $\sqrt{\sigma_{x}^{2}+\sigma_{y}^{2}}$ provides an excellent lower bound to the reported position errors as seen in Figure 17.

The degradation factors (ratio of reported error to expected error) in flux and position are correlated but not redundant, as seen in Figure 18.