gmm-recovery 0.1.2

GMM recovery impl by V.G.Spokoiniy and David

GMM-Recovery: Spokoiny's Method of Moments for Isotropic Gaussian Mixture Models

A lightweight Python implementation of the Generalized Method of Moments (GMM) for recovering the parameters of multivariate isotropic Gaussian mixture models, adhering strictly to the theoretical framework proposed by V. Spokoiny (2026).

1. Mathematical Formulation

Let $X_1, \dots, X_n \in \mathbb{R}^d$ be independent and identically distributed (i.i.d.) random vectors drawn from an isotropic Gaussian mixture with density:

$$f(x, \theta) = \sum_{k=1}^K \mu_k \sigma_k^{-d} \varphi\left(\frac{x - m_k}{\sigma_k}\right)$$

where $\varphi(x) = (2\pi)^{-d/2} \exp(-|x|^2 / 2)$, and the vector of unknown parameters is defined as:

$$\theta = {(m_k, \sigma_k, \mu_k), k = 1, \dots, K}, \quad \mu_k \geq 0, \quad \sum_{k=1}^K \mu_k = 1$$

1.1. System of Moments and Test Functions

Given a family of test functions ${\omega_j(x), j = 1, \dots, J}$, we define the empirical observables (moments):

$$Z_j = \frac{1}{n} \sum_{i=1}^n \omega_j(X_i), \quad j = 1, \dots, J$$

The test functions are constructed as isotropic Gaussian kernels centered at points $\xi_j \in \mathbb{R}^d$ with scale parameters $s_j > 0$:

$$\omega_j(x) = \mathcal{K}\left(\frac{x - \xi_j}{s_j}\right) = \exp\left( -\frac{|x - \xi_j|^2}{2s_j^2} \right)$$

The theoretical expectation $M_j(\theta) = \mathbb{E}[\omega_j(X)] = \int \omega_j(x)f(x, \theta)dx$ under the isotropic kernel yields a closed-form analytical expression (Lemma 1.1):

$$\psi_j(m_k, \sigma_k) = \left( 1 + \frac{\sigma_k^2}{s_j^2} \right)^{-d/2} \exp\left{ - \frac{|m_k - \xi_j|^2}{2(\sigma_k^2 + s_j^2)} \right}$$

$$M_j(\theta) = \sum_{k=1}^K \mu_k \psi_j(m_k, \sigma_k)$$

---

2. Alternating Minimization Procedure

To ensure numerical stability and eliminate the direct non-linear dependency on the moment residuals, we introduce an auxiliary image parameter $z \in \mathbb{R}^J$ and rewrite the extended objective function as (Equation 1.7):

$$\mathcal{L}(\theta, z) = \frac{1}{2} |Z - z|^2 + \frac{1}{2} |z - M(\theta)|^2$$

The optimization is solved via alternating minimization at each outer iteration $t$:

1. Step over $z$ (analytical): Updates the surrogate target by taking the average of empirical and estimated moments:

   $$z^{(t+1)} = \frac{1}{2}\left( Z + M(\theta^{(t)}) \right)$$

2. Step over $\theta$ (numerical): Minimizes the smooth quadratic discrepancy:

   $$\theta^{(t+1)} = \arg\min_{\theta} \frac{1}{2} |z^{(t+1)} - M(\theta)|^2$$

   This step is computed using the L-BFGS-B quasi-Newton optimizer with analytical gradients.

2.1. Parameter Constraints and Reparameterization

To enforce strict boundary conditions during the optimization steps:

• Mixture weights $\mu$ are parameterized via a softmax mapping: $\mu_k = e^{\eta_k} / \sum_l e^{\eta_l}$ for unconstrained parameters $\eta \in \mathbb{R}^K$.
• Standard deviations $\sigma$ are kept positive using an exponential map: $\sigma_k = e^{\lambda_k}$ for $\lambda_k \in \mathbb{R}$.

---

3. Test Function Selection Strategy

The library provides two strategies for generating the centers ${\xi_j}$ and scales ${s_j}$ of the test functions:

1. Data-dependent (Default): Test centers are perturbed versions of the actual data points, $\xi_j = X_j + \varepsilon_j$, where $\varepsilon_j \sim \mathcal{N}(0, \sigma_0^2 I_d)$. The bandwidths $s_j$ are adaptively set proportional to the median distance to the nearest neighbor. This concentrates moments in regions with high data density.
2. Multiscale (Optional): If a specific number of test functions $q$ is requested, they are sampled over a subset of the dataset and distributed across three bandwidth scales $s_j \in {0.3s_0, s_0, 3s_0}$ to capture both fine and coarse properties of the mixture structure.

---

4. Quick Start

Installation

pip install gmm_recovery