Prior Harnessing Index (PHI) for Galaxy Decomposition

• Prior Harnessing Index (PHI) is a Bayesian MCMC algorithm designed for the structural decomposition of galaxy images, ensuring physically plausible fits through adaptive prior management.
• It employs a triple-layer adaptive MCMC sampling strategy that refines proposal covariances and acceptance rates for efficient exploration of complex, multimodal parameter spaces.
• By rigorously constraining model parameters and excluding unphysical configurations, PHI enhances model selection and convergence diagnostics using metrics like BIC and Gelman–Rubin statistics.

The Prior Harnessing Index (PHI) is a fully Bayesian Markov-chain Monte Carlo (MCMC) algorithm developed for the structural decomposition of galaxy images. It combines informed prior selection with a triple-layer adaptive MCMC scheme to address the complex, multimodal parameter spaces inherent in fitting one- or two-component (bulge and disc) galaxy light profiles. PHI’s central innovation lies in the rigorous treatment of priors—termed “prior harnessing"—which prevent nonphysical model configurations, and in its robust exploration of the posterior distribution of galaxy structural parameters (Argyle et al., 2018).

1. Bayesian Construction and Likelihood Evaluation

PHI operates within a strict Bayesian framework to infer the parameter vector $\theta$—representing all structural and geometric parameters for one- or two-component models—given observed image data $D$, comprised of pixel values with known uncertainties $\sigma_i$. The likelihood for each pixel under the Gaussian noise assumption is

$$p(d_i \mid \theta) = (2\pi\,\sigma_i^2)^{-1/2}\,\exp\left\{-\frac{[d_i - f(x_i; \theta)]^2}{2\sigma_i^2}\right\},$$

yielding, for an image of $N$ pixels,

$$L(D \mid \theta) = \prod_{i=1}^N p(d_i \mid \theta).$$

The log-likelihood simplifies to

$$-2\ln L(D \mid \theta) = \sum_{i=1}^N \left[\frac{(d_i - f_i(\theta))^2}{\sigma_i^2} + \ln(2\pi \sigma_i^2)\right].$$

Priors adopted in the PHI framework are arbitrary and modular, enabling both individual and composite (joint) constraints. The Bayesian posterior is then proportional to $L(D|\theta)p(\theta)$, with priors expressed as

$$p(\theta) = \prod_j p_j(\theta_j) \times p_{\text{combo}}(\theta),$$

where $p_j$ denotes independent one-dimensional priors, and $D$0 encodes nonphysical exclusions, especially for two-component model reversals.

2. Surface Brightness Profile Parameterization

PHI models galaxy images using classical surface brightness laws:

• Sérsic Bulge: $D$1, where $D$2 is the Sérsic index and $D$3 is defined implicitly via $D$4 to ensure $D$5 encloses half the bulge luminosity.
• Exponential Disc: $D$6, with $D$7 the disc’s central intensity and $D$8 the scale length.

These prescriptions allow parameter inference for both single- and two-component model families, providing analytic tractability for the generative models fitted by the MCMC sampler.

3. Triple-Layer Adaptive MCMC Sampling

The central computational tool of PHI is its three-stage adaptive Metropolis MCMC procedure:

• Level 1: Blocked Adaptive Metropolis within Gibbs. Each parameter $D$9 is updated in turn by proposing $\sigma_i$0, where $\sigma_i$1 is the current parameter value. The step size $\sigma_i$2 is adapted every $\sigma_i$3 to target optimal acceptance rates ($\sigma_i$4 for single-parameter updates, $\sigma_i$5 for multidimensional blocks). Once all $\sigma_i$6 stabilize within $\sigma_i$7, adaptation ceases.
• Level 2: Adaptive Metropolis with Full Covariance. The sample covariance $\sigma_i$8 of preceding Level 1 draws is used to propose $\sigma_i$9, $$p(d_i \mid \theta) = (2\pi\,\sigma_i^2)^{-1/2}\,\exp\left\{-\frac{[d_i - f(x_i; \theta)]^2}{2\sigma_i^2}\right\},$$0 with $$p(d_i \mid \theta) = (2\pi\,\sigma_i^2)^{-1/2}\,\exp\left\{-\frac{[d_i - f(x_i; \theta)]^2}{2\sigma_i^2}\right\},$$1 the parameter dimension. Updates to $$p(d_i \mid \theta) = (2\pi\,\sigma_i^2)^{-1/2}\,\exp\left\{-\frac{[d_i - f(x_i; \theta)]^2}{2\sigma_i^2}\right\},$$2 continue every $$p(d_i \mid \theta) = (2\pi\,\sigma_i^2)^{-1/2}\,\exp\left\{-\frac{[d_i - f(x_i; \theta)]^2}{2\sigma_i^2}\right\},$$3 iterations, then are frozen when changes diminish.
• Level 3: Metropolis with Fixed Covariance. After fixing $$p(d_i \mid \theta) = (2\pi\,\sigma_i^2)^{-1/2}\,\exp\left\{-\frac{[d_i - f(x_i; \theta)]^2}{2\sigma_i^2}\right\},$$4 and $$p(d_i \mid \theta) = (2\pi\,\sigma_i^2)^{-1/2}\,\exp\left\{-\frac{[d_i - f(x_i; \theta)]^2}{2\sigma_i^2}\right\},$$5 from Level 2, multiple chains are run with proposals as in Level 2. Convergence diagnostics—Gelman–Rubin $$p(d_i \mid \theta) = (2\pi\,\sigma_i^2)^{-1/2}\,\exp\left\{-\frac{[d_i - f(x_i; \theta)]^2}{2\sigma_i^2}\right\},$$6 or Geweke’s Z-score—guide chain merging for the final posterior.

This hierarchical adaptation ensures efficient exploration of complex, correlated, and multimodal posteriors.

4. Prior Harnessing: Physical Constraints and Exclusion of Pathological Solutions

PHI’s approach to prior selection serves two main functions: (i) encoding minimal but essential physical knowledge, and (ii) explicitly excluding unphysical solutions, particularly for bulge-disc reversals in two-component fits.

• Parameter Priors: Key parameters—including $$p(d_i \mid \theta) = (2\pi\,\sigma_i^2)^{-1/2}\,\exp\left\{-\frac{[d_i - f(x_i; \theta)]^2}{2\sigma_i^2}\right\},$$7, $$p(d_i \mid \theta) = (2\pi\,\sigma_i^2)^{-1/2}\,\exp\left\{-\frac{[d_i - f(x_i; \theta)]^2}{2\sigma_i^2}\right\},$$8, $$p(d_i \mid \theta) = (2\pi\,\sigma_i^2)^{-1/2}\,\exp\left\{-\frac{[d_i - f(x_i; \theta)]^2}{2\sigma_i^2}\right\},$$9, $N$0—are taken as uniform in log-space over $N$1. The Sérsic index $N$2 has $N$3; axis ratio $N$4, $N$5; position angle $N$6; centroid $N$7 uniform over the image.
• Composite Priors: Imposed for two-component models:
  • Bulge/disc size ratio $N$8
  • Bulge-to-total flux $N$9
  • Bulge dominates inside $$L(D \mid \theta) = \prod_{i=1}^N p(d_i \mid \theta).$$0: $$L(D \mid \theta) = \prod_{i=1}^N p(d_i \mid \theta).$$1
  • Enforcing exactly one crossing of the projected 1D bulge and disc profiles, implemented via a Newton–Raphson root search.

These constraints robustly exclude parameter regions yielding swapped "inner" and "outer" model configurations, directly addressing pitfalls in classical two-component decomposition.

5. Diagnostics, Model Selection, and Goodness-of-Fit

PHI incorporates established Bayesian and frequentist diagnostics to ensure both fit robustness and correct model selection:

• Convergence: Level 3 employs the Gelman–Rubin $$L(D \mid \theta) = \prod_{i=1}^N p(d_i \mid \theta).$$2 statistic, where $$L(D \mid \theta) = \prod_{i=1}^N p(d_i \mid \theta).$$3 signals proper convergence. The Geweke diagnostic is optionally available for single chains.
• Model Selection via BIC: The Bayesian Information Criterion (BIC) distinguishes between one- and two-component fits. For models with $$L(D \mid \theta) = \prod_{i=1}^N p(d_i \mid \theta).$$4 and $$L(D \mid \theta) = \prod_{i=1}^N p(d_i \mid \theta).$$5 parameters over $$L(D \mid \theta) = \prod_{i=1}^N p(d_i \mid \theta).$$6 pixels:

$$L(D \mid \theta) = \prod_{i=1}^N p(d_i \mid \theta).$$7

Model preference is then quantified by $$L(D \mid \theta) = \prod_{i=1}^N p(d_i \mid \theta).$$8. In synthetic tests, $$L(D \mid \theta) = \prod_{i=1}^N p(d_i \mid \theta).$$9 yielded $$-2\ln L(D \mid \theta) = \sum_{i=1}^N \left[\frac{(d_i - f_i(\theta))^2}{\sigma_i^2} + \ln(2\pi \sigma_i^2)\right].$$0 completeness and $$-2\ln L(D \mid \theta) = \sum_{i=1}^N \left[\frac{(d_i - f_i(\theta))^2}{\sigma_i^2} + \ln(2\pi \sigma_i^2)\right].$$1 purity in separating single from bulge+disc galaxies. For real SDSS data, the separation is less distinct, necessitating possible supplementation with morphological priors or visual inspection.

6. Empirical Limitations and Reliability Considerations

PHI’s ability to recover unbiased, well-constrained parameters depends on specific galaxy and dataset properties:

• Bulge parameters exhibit substantially larger errors than disc parameters, especially when $$-2\ln L(D \mid \theta) = \sum_{i=1}^N \left[\frac{(d_i - f_i(\theta))^2}{\sigma_i^2} + \ln(2\pi \sigma_i^2)\right].$$2 (subdominant bulge flux) or bulges are marginally resolved with respect to the PSF (i.e., $$-2\ln L(D \mid \theta) = \sum_{i=1}^N \left[\frac{(d_i - f_i(\theta))^2}{\sigma_i^2} + \ln(2\pi \sigma_i^2)\right].$$3 PSF FWHM).
• High Sérsic index bulges ($$-2\ln L(D \mid \theta) = \sum_{i=1}^N \left[\frac{(d_i - f_i(\theta))^2}{\sigma_i^2} + \ln(2\pi \sigma_i^2)\right].$$4) show pronounced degeneracies in $$-2\ln L(D \mid \theta) = \sum_{i=1}^N \left[\frac{(d_i - f_i(\theta))^2}{\sigma_i^2} + \ln(2\pi \sigma_i^2)\right].$$5, causing bias in point estimates and slow convergence.
• Posterior covariances, particularly among $$-2\ln L(D \mid \theta) = \sum_{i=1}^N \left[\frac{(d_i - f_i(\theta))^2}{\sigma_i^2} + \ln(2\pi \sigma_i^2)\right].$$6, $$-2\ln L(D \mid \theta) = \sum_{i=1}^N \left[\frac{(d_i - f_i(\theta))^2}{\sigma_i^2} + \ln(2\pi \sigma_i^2)\right].$$7, and $$-2\ln L(D \mid \theta) = \sum_{i=1}^N \left[\frac{(d_i - f_i(\theta))^2}{\sigma_i^2} + \ln(2\pi \sigma_i^2)\right].$$8, are substantial and should be propagated in any downstream analysis.
• A plausible implication is that when bulge flux is weak or the PSF dominates, bulge parameter inference is systematically less robust, with potential impacts for population-level scaling relation studies.

7. Significance and Context within Galaxy Image Analysis

PHI formalizes and automates Bayesian photometric decomposition of galaxies, with key strengths arising from its explicit “prior harnessing” and adaptive posterior sampling. Its agreement with alternative codes on SDSS data, when controlling for input models and observational characteristics, demonstrates the validity of its Bayesian framework. The method’s capacity to avoid unphysical solutions, propagate parameter covariances, and offer statistically rigorous model differentiation represents significant progress over traditional optimization-based approaches in astrophysical image analysis (Argyle et al., 2018).