Multi-Gaussian Expansion (MGE)

• Multi-Gaussian Expansion (MGE) is a parametric technique that decomposes two-dimensional distributions into sums of elliptical or bivariate Gaussians for flexible imaging analysis.
• The method enables analytic operations like convolution, projection, and deprojection, which facilitates precise modeling in astrophysics including galaxy photometry and strong-lensing mass reconstruction.
• MGE applications extend to dynamical modeling via Jeans equations and even to image compression and face modeling, demonstrating versatility across various scientific domains.

The Multi-Gaussian Expansion (MGE) is a parametric decomposition technique for representing two-dimensional distributions as a sum of elliptical or general bivariate Gaussians. This framework allows for analytic manipulations—including convolution, projection, and deprojection—making MGE a highly flexible and robust model for a variety of scientific imaging and inference applications. First established in galaxy photometry and dynamics, MGE has seen successful deployment in galaxy surface-brightness decomposition, strong-lensing mass modeling, face image compression, and dynamical modeling via Jeans equations.

1. Mathematical Formalism and Parameterization

The core concept of MGE is the representation of an arbitrary two-dimensional function (e.g., surface brightness, mass density) as a linear sum of $N$ Gaussian components. Each component is parameterized by amplitude, spatial position/centroid, dispersion(s), and geometric orientation:

$$I(x, y) = \sum_{j=1}^{N} I_{0,j} \exp\left\{ -\frac{1}{2} \left[ \frac{x_j'^2}{\sigma_{x,j}^2} + \frac{y_j'^2}{\sigma_{y,j}^2} \right] \right\}$$

with

$$\begin{aligned} x_j' &= (x-x_j^c)\cos\theta_j + (y-y_j^c)\sin\theta_j \ y_j' &= -(x-x_j^c)\sin\theta_j + (y-y_j^c)\cos\theta_j \end{aligned}$$

Alternatively, for elliptical Gaussians with axis ratio $q_j = \sigma_{y,j}/\sigma_{x,j}$ and a single major-axis scale $\sigma_j$, the expression simplifies to

$$I(x, y) = \sum_{j=1}^{N} I_{0,j} \exp\left\{ -\frac{1}{2\sigma_{j}^2} \left[ x_j'^2 + \frac{y_j'^2}{q_j^2} \right] \right\}$$

The normalization can be chosen such that the integrated flux of each Gaussian is $I_{j}$, i.e., including a prefactor $1/(2\pi q_j \sigma_{j}^2)$ (Miller et al., 2021, He et al., 2024).

The functional form's linearity in $I_{0,j}$ and analytic properties under rotation, translation, and affine transforms enable powerful modeling and inference capabilities.

2. Analytic Operations: Convolution and Deprojection

One of MGE's main strengths is its preservation under analytic operations:

Convolution:

If the PSF is itself expressible as a sum of circular or elliptical Gaussians, convolution of the MGE model with the PSF is analytic: each model–PSF Gaussian pair yields another Gaussian, with summed variances and axis ratios transformed by the PSF's geometry. For $M$ PSF components,

$$I(x, y) = \sum_{j=1}^{N} I_{0,j} \exp\left\{ -\frac{1}{2} \left[ \frac{x_j'^2}{\sigma_{x,j}^2} + \frac{y_j'^2}{\sigma_{y,j}^2} \right] \right\}$$0

where $$I(x, y) = \sum_{j=1}^{N} I_{0,j} \exp\left\{ -\frac{1}{2} \left[ \frac{x_j'^2}{\sigma_{x,j}^2} + \frac{y_j'^2}{\sigma_{y,j}^2} \right] \right\}$$1 and $$I(x, y) = \sum_{j=1}^{N} I_{0,j} \exp\left\{ -\frac{1}{2} \left[ \frac{x_j'^2}{\sigma_{x,j}^2} + \frac{y_j'^2}{\sigma_{y,j}^2} \right] \right\}$$2 (Miller et al., 2021, He et al., 2024).

Projection/Deprojection:

Given an inclination $$I(x, y) = \sum_{j=1}^{N} I_{0,j} \exp\left\{ -\frac{1}{2} \left[ \frac{x_j'^2}{\sigma_{x,j}^2} + \frac{y_j'^2}{\sigma_{y,j}^2} \right] \right\}$$3 and symmetry assumptions (e.g., oblate axisymmetry), each 2D Gaussian component uniquely deprojects to a 3D Gaussian. The intrinsic axis ratio is

$$I(x, y) = \sum_{j=1}^{N} I_{0,j} \exp\left\{ -\frac{1}{2} \left[ \frac{x_j'^2}{\sigma_{x,j}^2} + \frac{y_j'^2}{\sigma_{y,j}^2} \right] \right\}$$4

and the 3D luminosity/mass density

$$I(x, y) = \sum_{j=1}^{N} I_{0,j} \exp\left\{ -\frac{1}{2} \left[ \frac{x_j'^2}{\sigma_{x,j}^2} + \frac{y_j'^2}{\sigma_{y,j}^2} \right] \right\}$$5

enabling full dynamical modeling and potential theory (Li et al., 2015).

3. Methodologies for Fitting and Inference

Bayesian Image Modeling (Miller et al., 2021)

MGE models are fitted to imaging data by maximum-likelihood or full posterior-sampling techniques. Typically:

• Global structural parameters $$I(x, y) = \sum_{j=1}^{N} I_{0,j} \exp\left\{ -\frac{1}{2} \left[ \frac{x_j'^2}{\sigma_{x,j}^2} + \frac{y_j'^2}{\sigma_{y,j}^2} \right] \right\}$$6 are optimized via least-squares (e.g., Trust-Region Reflective, Levenberg–Marquardt).
• With positions and geometry fixed, the amplitudes ($$I(x, y) = \sum_{j=1}^{N} I_{0,j} \exp\left\{ -\frac{1}{2} \left[ \frac{x_j'^2}{\sigma_{x,j}^2} + \frac{y_j'^2}{\sigma_{y,j}^2} \right] \right\}$$7) and sky parameters become the main parameters for Bayesian inference.
• Empirical Gaussian priors, centered on best-fit values, facilitate rapid convergence during nested sampling or MCMC, with dynamic nested sampling (e.g., dynesty) commonly employed.
• For efficiency, the PSF-convolved images of each unit-flux component $$I(x, y) = \sum_{j=1}^{N} I_{0,j} \exp\left\{ -\frac{1}{2} \left[ \frac{x_j'^2}{\sigma_{x,j}^2} + \frac{y_j'^2}{\sigma_{y,j}^2} \right] \right\}$$8 are pre-rendered, and the model at any likelihood evaluation is $$I(x, y) = \sum_{j=1}^{N} I_{0,j} \exp\left\{ -\frac{1}{2} \left[ \frac{x_j'^2}{\sigma_{x,j}^2} + \frac{y_j'^2}{\sigma_{y,j}^2} \right] \right\}$$9, reducing the computational burden by $$\begin{aligned} x_j' &= (x-x_j^c)\cos\theta_j + (y-y_j^c)\sin\theta_j \ y_j' &= -(x-x_j^c)\sin\theta_j + (y-y_j^c)\cos\theta_j \end{aligned}$$0.

Semi-Linear Inversion with Lens Modeling (He et al., 2024)

For strong-lensing images, MGE lens-light modeling is integrated with semi-linear inversion for pixelized source reconstruction. The data model is

$$\begin{aligned} x_j' &= (x-x_j^c)\cos\theta_j + (y-y_j^c)\sin\theta_j \ y_j' &= -(x-x_j^c)\sin\theta_j + (y-y_j^c)\cos\theta_j \end{aligned}$$1

where $$\begin{aligned} x_j' &= (x-x_j^c)\cos\theta_j + (y-y_j^c)\sin\theta_j \ y_j' &= -(x-x_j^c)\sin\theta_j + (y-y_j^c)\cos\theta_j \end{aligned}$$2 accounts for lens deflections, PSF, etc. Amplitudes of both the source pixels and MGE Gaussians are solved simultaneously, regularized, and constrained via non-negativity to avoid overfitting lens light into the lensed source image. Bayesian evidence is maximized for model selection and exploration.

Other Applications: GmFace (Zhang et al., 2020)

For modeling facial images, the "GmFace" framework expresses intensity as

$$\begin{aligned} x_j' &= (x-x_j^c)\cos\theta_j + (y-y_j^c)\sin\theta_j \ y_j' &= -(x-x_j^c)\sin\theta_j + (y-y_j^c)\cos\theta_j \end{aligned}$$3

with all parameters learned via gradient descent in a GmNet architecture, directly mapping network weights to model parameters.

4. Validation, Performance, and Systematics

Extensive validation is reported across domains:

• Galaxy profile inference: MGE models recover total fluxes to $$\begin{aligned} x_j' &= (x-x_j^c)\cos\theta_j + (y-y_j^c)\sin\theta_j \ y_j' &= -(x-x_j^c)\sin\theta_j + (y-y_j^c)\cos\theta_j \end{aligned}$$4 mag and effective radii to $$\begin{aligned} x_j' &= (x-x_j^c)\cos\theta_j + (y-y_j^c)\sin\theta_j \ y_j' &= -(x-x_j^c)\sin\theta_j + (y-y_j^c)\cos\theta_j \end{aligned}$$5 for injected mock galaxies in HST images (Miller et al., 2021).
• Lens modeling: Typical residuals in lens-light subtraction remain below $$\begin{aligned} x_j' &= (x-x_j^c)\cos\theta_j + (y-y_j^c)\sin\theta_j \ y_j' &= -(x-x_j^c)\sin\theta_j + (y-y_j^c)\cos\theta_j \end{aligned}$$6 everywhere except Poisson-dominated centers, and fit residuals consistently lie within $$\begin{aligned} x_j' &= (x-x_j^c)\cos\theta_j + (y-y_j^c)\sin\theta_j \ y_j' &= -(x-x_j^c)\sin\theta_j + (y-y_j^c)\cos\theta_j \end{aligned}$$7 of the noise (He et al., 2024).
• Dynamical modeling: Using the JAM (Jeans Anisotropic MGE) technique, total mass within $$\begin{aligned} x_j' &= (x-x_j^c)\cos\theta_j + (y-y_j^c)\sin\theta_j \ y_j' &= -(x-x_j^c)\sin\theta_j + (y-y_j^c)\cos\theta_j \end{aligned}$$8 is recovered to $$\begin{aligned} x_j' &= (x-x_j^c)\cos\theta_j + (y-y_j^c)\sin\theta_j \ y_j' &= -(x-x_j^c)\sin\theta_j + (y-y_j^c)\cos\theta_j \end{aligned}$$9 scatter, with $q_j = \sigma_{y,j}/\sigma_{x,j}$0 bias $q_j = \sigma_{y,j}/\sigma_{x,j}$1 for oblate galaxies; higher scatter (and bias) is seen for prolate/triaxial systems (Li et al., 2015).

Identified Limitations and Mitigations:

• Insufficient number of Gaussians can introduce oscillatory "wiggles" or misfit outer/inner slopes. Bayesian model averaging across several component sets mitigates biases from discretization (Miller et al., 2021).
• Fixed $q_j = \sigma_{y,j}/\sigma_{x,j}$2 per Gaussian may not capture radially varying isophotes (e.g., twisting), requiring multiple sets with independently varying geometry (He et al., 2024).
• In dynamical applications, degeneracy between $q_j = \sigma_{y,j}/\sigma_{x,j}$3 and $q_j = \sigma_{y,j}/\sigma_{x,j}$4 is pronounced, especially at low resolution, requiring very high-fidelity imaging and careful regularization.
• Over-fitting is possible unless non-negativity of amplitudes is enforced, particularly in blended or strongly lensed systems (He et al., 2024).

5. Principal Domains of Application

Application Domain	MGE Role	Benchmark Reference
Galaxy photometry	Surface-brightness decomposition, PSF convolution	(Miller et al., 2021, Li et al., 2015)
Strong gravitational lensing	Lens light & mass modeling, isophotal structure	(He et al., 2024)
Stellar dynamics/JAM	Axisymmetric/oblate mass modeling, potential theory	(Li et al., 2015)
Face image modeling (GmFace)	Image compression, transformable representation	(Zhang et al., 2020)

In galactic contexts, MGE enables rapid, high-fidelity modeling of light distributions without reliance on fixed parametric forms (e.g., Sérsic), analytic treatment of PSF, and straightforward deprojection for dynamical modeling. In strong lensing, it facilitates the separation of lens and source light and uncovers subtle isophotal twist or boxiness. In image-processing contexts, parameter manipulations allow geometric transforms (scaling, rotation, translation) directly in parameter space.

6. Implementation Recommendations and Limitations

• Number of Gaussians is typically $q_j = \sigma_{y,j}/\sigma_{x,j}$5–$q_j = \sigma_{y,j}/\sigma_{x,j}$6, logarithmically spaced in dispersion to sample inner to outer galaxy regions. Fewer components lead to systematic features, more add computation with marginal benefit (Miller et al., 2021, Li et al., 2015).
• Robust least-squares optimization of geometric parameters followed by Bayesian or semi-linear fitting of amplitudes is standard (Miller et al., 2021, He et al., 2024).
• Regularization (smoothness, monotonicity) on amplitude vectors suppresses overfitting noise (He et al., 2024, Li et al., 2015).
• In dynamical models, highest-available spatial resolution is required to avoid central mass–light degeneracy; restrict analysis to oblate/mildly triaxial systems for unbiased recovery (Li et al., 2015).
• Non-negativity constraints are critical in blended-source contexts to prevent unphysical solutions (He et al., 2024).
• Discretization and parameter correlations (especially for neighboring dispersion scales) can impact inference; Bayesian model averaging addresses these for uncertainty quantification (Miller et al., 2021).
• Empirical priors, Hessian-derived uncertainties, and log-uniform alternatives influence convergence and uncertainty estimates.

7. Extensions and Future Directions

Recent formulations extend MGE to support large-scale joint inference (e.g., lens mass + source structure), multicomponent galaxies (stellar plus dark matter), and generalized isophote modeling (boxiness, disky features, radial twists) (He et al., 2024). Open-source tools such as imcascade (Miller et al., 2021) and public codes in the lensing community (e.g., PyAutoLens) natively implement the MGE formalism.

A plausible implication is wider adoption in multiwavelength and multimodal imaging domains, where analytic transformability and efficient posterior sampling are essential. The capacity for transparent geometric transformation of components (as in GmFace), and analytic manipulation for convolution and projection, position MGE as a foundational tool for quantitative image analysis across astrophysics, computer vision, and data-driven physical modeling.