Mixtures Closest To A Given Measure: A Semidefinite Programming Approach

Summary

The paper studies a classical unsupervised problem — approximating an unknown probability measure by a mixture from a parametric family (Gaussian, Poisson, exponential) — under the constraint that the target is accessed only through its moments and that mixture parameters lie in a compact basic semi-algebraic set. Building on Lasserre (2024, 2025), the authors formulate the $W_2$- and TV-distance approximation problems as generalised moment problems over a mixing measure $\phi$, add a polynomial regulariser $\varepsilon L_\phi(R)$, and construct Moment-SOS SDP hierarchies whose relaxations converge asymptotically to the optimal value; a flat-extension rank condition yields finite convergence and an atomic mixing measure from which the mixture order $\widehat K$ and component parameters can be extracted by GNS-style multiplication matrices. The conceptual move — treating the mixing measure itself as the decision variable and reading $\widehat K$ off the rank of an SDP moment matrix — is elegant and connects the mixture-order and parameter-estimation problems in a single convex-relaxation pipeline, sidestepping the non-convexity that plagues EM/MLE and the ad-hoc order-selection heuristics that dominate the literature. Lemma A.1, a multivariate Carleman bound for Gaussian, Poisson and exponential mixtures, is a genuinely useful side result. On the applied side, the authors demonstrate the pipeline as an initialiser that reduces the number of $k$-means/EM iterations on synthetic Gaussian mixtures ($K\in\{2,5\}$ in $\mathbb R^2$) and on MNIST after PCA-to-2D. Independent re-execution confirms the univariate uniform illustration at the aggregate level (both $\widehat K=4$ and $\tau\approx 0$ reproduce numerically), but the recovered component weights and most of the multivariate empirical results could not be verified from the shipped code. The main conceptual and methodological limitation is empirical scope and scalability: all reported experiments are in $n\le 2$ ambient dimensions, the SDP size grows quickly with $n$ and $p$, the compactness of $S_\theta$ (called "non-restrictive") is in fact a substantive modelling assumption, and the univariate route proposed for higher dimensions ends in an acknowledged combinatorial explosion without a resolution. The theoretical framework is solid and the SDP formulation is a real conceptual advance, but the paper overstates originality and optimality in places, and the released code only partially supports the applied claims.

Strengths

• Conceptual contribution. Formulating mixture approximation as a generalised moment problem in the mixing measure and lifting it to a Moment-SOS SDP hierarchy is a clean, principled way to jointly recover the mixture order and the component parameters, and it plausibly generalises well beyond the illustrations shown.
• Convergence guarantees. Theorems 3.1, 3.2, A.1 and A.2 establish existence of relaxation optima, asymptotic convergence, and finite convergence under the flat-extension condition, with sensible extraction (Algorithm 2) via multiplication matrices — a rigorous package for a problem where guarantees are usually weak.
• Independent lemma. The multivariate Carleman-type Lemma A.1 for Gaussian, Poisson and exponential mixtures is a clean auxiliary result that other moment-based mixture work can reuse.
• TV via Hahn–Jordan. Turning TV minimisation into a moment-domination LP via the Hahn–Jordan decomposition and translating the domination $\psi_-\le\nu_\phi$ into a moment-matrix inequality is a nice and non-obvious tractability device.
• Practical initialiser. The $K=2$ and $K=5$ experiments and the MNIST panel (a) show that SDP-extracted means yield meaningfully fewer $k$-means/EM iterations than random starts, supporting the "preprocessing" positioning; independent re-execution corroborates the univariate SDP core, with $\widehat K = 4$ and a numerically zero relaxation optimum on the uniform target.
• Honest limitations. The conclusion openly acknowledges scalability as the main bottleneck and points to concrete future directions (symmetry exploitation, sliced formulations).

Weaknesses

• "Non-restrictive" compactness. Every convergence proof and the Gaussian example itself rely on compactness of $S_\theta$ (bounded means, $0<a\le\Sigma\le b\,\mathbf I$). For unbounded Gaussians/Poissons/exponentials this is a real modelling choice — how one chooses the enclosing ball radius $R$ without prior knowledge of the mixture support is not discussed.
• Overstated originality. The claim that polynomial-in-parameter moments are "overlooked in the optimization and machine-learning literature" is at odds with the method-of-moments and tensor-decomposition papers the authors themselves cite (Améndola et al.; Wu & Yang; Khouja et al.; Pereira et al.; Lindberg et al.). The genuinely novel piece is the SDP lift over the mixing measure.
• Optimality claim vs. EM. The statement that EM-recovered mixtures "cannot be better" than Algorithm 1 in the uniform illustration presumes finite convergence at $d=4$ and equates optimality on $\mu$ with optimality on $\mu^N$ up to eight moments; this deserves qualification (target vs. empirical; flat-extension check for this instance).
• Empirical scope. The $K=5$ Gaussian experiment fixes a single $(c,\epsilon(\Sigma))=(5,0.25)$; the 90%/67% and 31.8%/36.5% headline improvements are anchored to one point in the space that separability and eccentricity were introduced to parameterise. MNIST is only tested after PCA-to-2D, and no experiment goes beyond $n=2$.
• Missing baselines. The $K=2$/$K=5$/MNIST experiments compare only against random (or, in the released MNIST notebook, $k$-means++) initialisation, not against BIC-, silhouette-, or gap-statistic-based order selection with an equally warm start. Since one selling point is that the SDP simultaneously delivers $\widehat K$, at least one such baseline is expected.
• Uniform illustration $\widehat K=4$. The illustration on $\mathrm{Unif}[0,1]$ returns exactly four components for both distances, matching the relaxation order $d=4$. Whether this is coincidence, a rank-estimator artefact of the moment-matrix size at $d=4$, or a genuine property of the target is not analysed.
• Table 1 order selection. Remark B.1's own mode rule gives $\widehat K=4$ at $\varepsilon=10^{-5}$ and $\widehat K=3$ at $\varepsilon=0.1$; the reported "recovery of ground truth $K=3$" then depends on choosing the matching $\varepsilon$. A principled $\varepsilon$-selection rule is missing, and the "weak $\varepsilon$ overfits, strong $\varepsilon$ recovers true $K$" mechanism is drawn from one dataset with two $\varepsilon$ values.
• Univariate stitching gap. Remark B.1 concedes that the univariate route yields combinatorially many candidate multivariate mixtures with no pruning strategy; since this is the main proposal for $n>2$, the practical utility of the pipeline in higher dimensions is left unresolved. This is compounded by Table 1 already showing $\widehat K_i$ spread across $\{1,\dots,5\}$.
• Scalability not experimentally probed. Section 5 identifies scalability as the main limitation, but no experiment sweeps $n$ or $p$; the reader gets solve-time information at $n=2$ only.
• TV domination with empirical moments. The block $\mathbf M_d(\mu)\succeq\mathbf M_d(\psi_+)$ requires the moment matrix of the true target, but only $\mathbf M_d(\mu^N)$ is available; the paper does not discuss how sampling noise interacts with this PSD domination or whether a PSD projection of $\mathbf M_d(\mu^N)$ is needed.
• Consistency of Algorithm 1. The pseudocode uses an exact rank equality for the flatness gate, then estimates $\widehat K$ with a cumulative-energy heuristic and warns that exact rank is numerically unreliable. The two notions of "rank" need to be reconciled.
• Technical issues in Appendix A. The two definitions of the truncation degrees $d_j$ (main body $\lceil\deg(r_j)/2\rceil$ vs. Appendix A.1 $\deg(2r_j)$) disagree; the quadratic module $Q(r)$ is written as a subset of $\mathbb R[\mathbf x]$ though the constraints live in $\mathbb R[\theta]$; Slater's condition for the TV pair is asserted without exhibiting a strictly feasible point; the domination $\psi_-\le\nu_\phi\Leftrightarrow \mathbf M_d(p;\phi)\succeq\mathbf M_d(\psi_-)$ is stated as a two-way equivalence when only the forward direction is generally valid; the support of a polynomial is defined with $f_\alpha=0$ where $f_\alpha\ne 0$ is intended.
• Illustration parameter set. $S_{m,\sigma}=\{(m,\sigma):m\ge m^2,\sigma\ge\sigma^2\}=[0,1]^2$ admits $\sigma=0$ (Dirac components), inconsistent with the strict positivity used for $\Sigma$ elsewhere and potentially relevant to why $\widehat K=4$ is recovered.
• PCA-to-2D on MNIST. The claim that PCA to two components "preserves the essential cluster structure" is asserted without a variance-retention or per-class-separation number; it also means the MNIST experiment does not stress-test the scalability limitation.

Reproducibility & code

Independent re-execution — with the SDP solver swapped from Mosek to the open-source Clarabel, since Mosek requires a commercial licence — reaches the univariate uniform experiment in Appendix C and cannot reach the multivariate and MNIST experiments. Findings:

• Headline aggregates on the uniform target reproduce. The regularised W2 and TV optima at $d=4$, $\varepsilon=0$ come out at $\approx -2.8\times 10^{-9}$ and $\approx -3.6\times 10^{-9}$, both numerically zero and matching the paper's $\tau^{W_2}_{4,0,R}=\tau^{TV}_{4,0,R}=0$. The cumulative-energy estimator returns $\widehat K = 4$ for both distances, matching Appendix C.
• Recovered component parameters partially reproduce. The four recovered means agree with the paper to within $\sim 0.02$ and the standard deviations to within $\sim 0.011$ (once a variance-vs-$\sigma$ storage convention is normalised), but the recovered weights deviate substantially (e.g. the broad W2 component's weight is $\approx 0.14$ against the paper's $0.212$). Appendix C reports these weights to eight significant digits, so this discrepancy is not just noise.
• Missing driver for the Figure 2 illustration. The Mixtures_Via_SDPs repository contains the univariate builders, extract_CF, and rank_by_energy, but no end-to-end notebook that draws the $N=2000$ Uniform$[0,1]$ samples, solves the two relaxations at $\varepsilon=0$, prints the objectives, extracts the four W2/TV components, and produces Figure 2 with the BIC/CV EM baselines. Reproducing the illustration required re-authoring the driver from scratch.
• Multivariate and Table 1 experiments not exercised end-to-end. The K=5, MNIST digits $\{0,1,2\}$, and Table 1 experiments were not reached in this reproduction: the shipped code has no Project.toml/Manifest.toml pinning TSSOS/DynamicPolynomials versions, has no open-source-solver path, and the K=5 aggregate percentages (90%/67%, 31.8%/36.5%, 18.63%/13.91%) are quoted in the text but not computed inside any released notebook.
• Timings not instrumented. No @time/@elapsed calls surround the SDP assembly, SDP solve, or extract_CF, so the 0.1755 s / 0.1506 s / 0.00036 s averages cannot be checked against the code without adding the instrumentation.

Recommended Changes

Essential

• Soften compactness and originality wording. Replace "non-restrictive" with a discussion of when the compact $S_\theta$ can be justified a priori (and how to choose $R$ otherwise), and rephrase the "overlooked" claim to point specifically at the SDP lift over the mixing measure (Weaknesses: "Non-restrictive" compactness; Overstated originality).
• Qualify the "cannot be better than Alg. 1" statement. State that the guarantee is with respect to the empirical target $\mu^N$ and the truncated moment problem, and either verify the flat-extension condition for the uniform illustration or explicitly note verification as pending (Weaknesses: Optimality claim vs. EM).
• Add a $(c,\epsilon(\Sigma))$ sweep and a model-selection baseline. Report iteration-reduction numbers on a grid of $(c,\epsilon(\Sigma))$ values for the $K=5$ experiment, and add a BIC- or gap-statistic-based order selector paired with $k$-means++/EM so the joint $(\widehat K,\text{init})$ contribution of the SDP is measurable (Weaknesses: Empirical scope; Missing baselines).
• Add a scalability probe. Include even a small experiment showing solve time and quality of $\widehat K$ as $n$ (or $p$) grows from 2 to, say, 10 on Gaussian mixtures, directly addressing the stated main limitation (Weaknesses: Scalability not experimentally probed).
• Unify the two rank criteria in Algorithm 1. Replace the exact-rank equality with a numerical-rank check tied to tol, matching the cumulative-energy $\widehat K$ estimator that is actually used (Weaknesses: Consistency of Algorithm 1).
• Fix the Appendix A technical issues. Reconcile $d_j$ with the main-body $\lceil\deg(r_j)/2\rceil$ convention, restate $Q(r)$ as a subset of $\mathbb R[\theta]$, exhibit a strictly feasible triple for the TV pair (or replace the Slater assertion with a proof sketch), weaken $\Leftrightarrow$ to $\Rightarrow$ in $\psi_-\le\nu_\phi\Rightarrow \mathbf M_d(p;\phi)\succeq\mathbf M_d(\psi_-)$ unless the additional determinacy hypothesis is stated, and fix the polynomial-support definition to use $f_\alpha\ne 0$ (Weaknesses: Technical issues in Appendix A).
• Ship the uniform-[0,1] illustration driver. Add a script that samples $N=2000$ from $\mathrm{Unif}[0,1]$, solves the $W_2$ and TV relaxations at $\varepsilon=0$, prints $\tau$, extracts the four components listed in Appendix C, and fits the BIC/CV EM baselines used in Figure 2 (Reproducibility & code: Missing driver for the Figure 2 illustration).
• Pin dependencies and provide an open-source-solver path. Add a Project.toml/Manifest.toml fixing the TSSOS and DynamicPolynomials versions actually used, and either document a Clarabel/SCS pipeline or explicitly state that Mosek is required (Reproducibility & code: Multivariate and Table 1 experiments not exercised end-to-end).
• Reconcile recovered weights. State which extraction routine and numerical settings produced the eight-significant-digit weight vectors in Appendix C, so third parties obtain the same $\phi$ rather than only the same $\widehat K$ and means (Reproducibility & code: Recovered component parameters partially reproduce).

Suggested

• Address the univariate stitching combinatorial explosion. Describe or evaluate a pruning strategy (e.g. a joint moment-consistency criterion, or the random-projection variant of Remark B.2) so the univariate route actually delivers a bounded set of multivariate candidates (Weaknesses: Univariate stitching gap).
• Sensitivity of $\widehat K$ to $d$ and tol on the uniform target. Report $\widehat K$ as a function of $d$ so readers can distinguish a genuine four-component fit from an artefact of the moment-matrix size (Weaknesses: Uniform illustration $\widehat K=4$).
• Extend the $\varepsilon$-sweep in Table 1. Report at least a fine $\varepsilon$-grid and, ideally, one additional dataset, before presenting "weak $\varepsilon$ overfits, strong $\varepsilon$ recovers true $K$" as a mechanism (Weaknesses: Table 1 order selection).
• Empirical-moment PSD safeguard. Either add code that projects $\mathbf M_d(\mu^N)$ onto the PSD cone before enforcing $\mathbf M_d(\mu^N)\succeq \mathbf M_d(\psi_+)$, or characterise the $N$ regime in which the projection is negligible (Weaknesses: TV domination with empirical moments).
• Aggregate the K=5 percentages inside the notebook and instrument timings. Add the lines that compute 90%/67%, 31.8%/36.5%, 18.63%/13.91% directly from the stored per-instance arrays, and place @elapsed around the SDP assembly/solve and extract_CF calls (Reproducibility & code: Multivariate and Table 1 experiments not exercised end-to-end; Timings not instrumented).
• Correct minor notation and set descriptions. Tighten the $S_{m,\sigma}$ set so $\sigma=0$ is excluded, and quantify PCA-to-2D's variance retention and per-class separation on MNIST (Weaknesses: Illustration parameter set; PCA-to-2D on MNIST).