Connecting Score Matching, Maximum Likelihood, and Expectation-Maximization in Mixed Linear Regression

Zhankun Luo; Abolfazl Hashemi

cs.LGSep 4, 2026

Connecting Score Matching, Maximum Likelihood, and Expectation-Maximization in Mixed Linear Regression

Authors: Zhankun Luo, Abolfazl Hashemi

Organizations: School of Electrical and Computer Engineering Purdue University, West Lafayette, IN, USA

Abstract

We study variance-preserving diffusion of the response in mixed linear regression (MLR) with unknown mixing weights. Our analysis separates the statistical guarantees of score matching from the loss geometry and optimization signal at a fixed diffusion noise level. The KL divergence links the denoising score matching objective integrated over the diffusion path with the likelihood and a terminal discrepancy. Under mild regularity conditions and terminal schedule, the resulting estimator converges up to the ground truth parameters of MLR, and its scaled error converges to the Gaussian limit of the maximum-likelihood estimator. At a fixed scale of the diffusion noise level, we derive a decomposition linking the score matching loss to cross-entropy and Expectation-Maximization (EM) operators. This decomposition yields an EM-related low-noise gradient expansion with additional correction terms of latent variance. In the high-noise limit, we further characterize gradient descent on this limiting loss under isotropic covariance. Along fixed high signal-to-noise ratio rays, the score matching imbalance gradient and the latent-variance term tend to zero pointwise. Numerical experiments illustrate our theoretical findings and statistical guarantees.

Figures & tables

Figure 1: Illustration of the finite-sample estimation results in Theorem 3.2 . The experiment uses a two-parameter Gaussian 2MLR instance, the theorem’s joint sign alignment, and a terminal sequence whose terminal-scale condition is verified in Appendix E . Panel (a) reports the median and 90% quantile of the error after joint sign alignment in the Fisher metric, together with 95% bootstrap intervals and the corresponding limiting quantiles of $\sqrt{\chi^{2}_{2}/n}$ . Over the three largest sample sizes, the fitted slopes are $-0.471$ and $-0.526$ , close to the predicted $-1/2$ in the rate of $n^{-1/2}$ . Panel (b) compares the squared, $n$ -scaled Fisher error with $\chi^{2}_{2}$ at the smallest and largest sample sizes. Coverage of the joint 95% Fisher ellipse ranges from 93.5% to 96.5%, equals 95.0% at the largest sample size, and has a corresponding covariance discrepancy of 0.070. These diagnostics provide finite-sample illustrations under the alignment and terminal-scale conditions of Main Theorem 3.2 . Appendix E provides additional calibration details.

Figure 2: First-order scaling of the endpoint gradient remainders in Propositions 4.7 and 4.8 . At one fixed model pair, the panels plot componentwise remainders for the complete low-noise expression in Proposition 4.7 and the effective high-noise coefficient expression in Proposition 4.8 . The fitted log–log slopes track the guides with unit slope, consistent with the respective

\mathcal{O}(1-\bar{\alpha}_{t})

and

\mathcal{O}(\bar{\alpha}_{t})

bounds. The panels report one selected instance, whose settings and fitted values appear in Appendix E .

Figure 3: Optimization Dynamics of Gradient Descent on the High-Noise Objective in Theorem 5.1 . The gradient descent update with a fixed step size is iterated from initializations satisfying the admissibility condition in Theorem D.1 , and every displayed path preserves that condition and has nonincreasing loss. (a) For the degenerate target, the computed loss remains below the $K^{-2}$ upper bound in Theorem 5.1 , and its fitted log–log tail slope is $-2-1.50\times 10^{-2}$ , tracking the anchored $K^{-2}$ guide over the displayed tail. (b) For $c^{\ast}\neq\mathbf{0}$ , the three displayed generic initializations lie outside the exceptional saddle basin and show initialization-dependent transients followed by geometric loss decay, consistent with the $\mathcal{O}(\log(1/\epsilon))$ conclusion. (c) The phase portrait in $(\langle\theta_{K},c^{\ast}\rangle/\|c^{\ast}\|,\nu_{K})$ coordinates uses open circles and filled diamonds to mark the initial and final iterates. The dotted curves are the two branches of the global-minimum set $\{(\theta,\nu):\theta\tanh\nu=c^{\ast}\}$ . By contrast, the saddle set $\{(\theta,0):\langle\theta,c^{\ast}\rangle=0\}$ projects to the single point $(0,0)$ in these coordinates and is marked by the black cross. The observed degenerate slope diagnoses finite-window agreement with the $K^{-2}$ upper bound in Theorem 5.1 . Appendix E provides the configurations for generating the figures.

Figure 4: Score Matching Blindness at High SNR in Proposition 5.2 . The displayed instance uses $d=1$ , $X\sim\mathcal{N}(0,1)$ , $\bar{\alpha}_{t}=\tfrac{1}{2}$ , the fixed candidate $(\theta,\nu)=(0.80,-0.35)$ , and $\nu^{\ast}=0.55$ , while $\theta^{\ast}>0$ grows along the positive ray. This setting satisfies Proposition 5.2 . (a)–(b) Following Wenliang and Kanagawa (2020) , samples from the varying truths are overlaid on the fixed-candidate sensitivity $|\partial_{\nu}s_{\theta_{t},\nu}(y_{t},x)|=|\mu_{t}|\sech^{2}(\mu_{t}y_{t}+\nu)$ on a log scale. The normalized responses concentrate on $y_{t}/\|\theta^{\ast}\|\simeq\pm\sqrt{\bar{\alpha}_{t}}\,x$ , and the appreciable sensitivity concentrates near the ray intersection. (c) Quadrature shows that $|\nabla_{\nu}\mathcal{L}_{t}|$ and the latent variance term $V_{t}=V(\theta_{t},\nu)/\bar{\alpha}_{t}$ decrease along the fixed ray, providing pointwise evidence at one fixed model and scale. Appendix Figure 6 provides fixed-scale and cross-entropy controls together with an auxiliary moving-candidate diagnostic for a separate setting. At high noise, Proposition 4.8 gives $\frac{1}{2}\|\theta\tanh\nu-c^{\ast}\|_{\Sigma}^{2}$ , whose directional signal vanishes at balance. At low noise, Proposition 4.7 retains operator-Jacobian and posterior-variance corrections beyond the EM residuals.

Appendix figures & tables9 assets

Supplementary material from the paper’s appendix.

Appendix

Quantity	Use in the experiments
$(\theta,\nu)$ and $(\theta^{\ast},\nu^{\ast})$	Model and ground-truth parameters in every study.
$\bar{\alpha}_{t}$ and $\bar{\alpha}_{T_{n}}$	Scale of the diffusion noise level in the endpoint and blindness studies, and terminal scale in the finite-sample study.
$(\widehat{\theta}_{n,T_{n}}^{\,\mathsf{SM}},\widehat{\nu}_{n,T_{n}}^{\,\mathsf{SM}})$ and $s_{n}$	Integrated score matching estimator and the joint sign alignment in Theorem 3.2 .
$I(\theta^{\ast},\nu^{\ast})$	Fisher-information benchmark for the standardized aligned estimator.
$M$ , $N$ , $V$ , and $V_{t}$	Terms in the complete low-noise reference and the blindness study.
$\Sigma$	Covariate covariance in the endpoint and limiting-loss calculations.

Appendix

Table 1: Experimental notation. Main-paper quantities reused in the numerical experiments.

Study	Model and numerical setting
Finite sample	$d=1$ , $(\theta^{\ast},\nu^{\ast})=(1.10,0.55)$ , $n=250\cdot 2^{j}$ for $j=0,\ldots,4$ , with 200 independent repetitions at each $n$ , $\bar{\alpha}_{0}=1$ , and $\bar{\alpha}_{T_{n}}=n^{-3/2}$ . The diffusion integral uses 24-point Gauss–Legendre quadrature, the expectation over $\xi$ uses 11-point Gauss–Hermite quadrature, and each fit uses four starts. The Fisher information uses an 80-point Gauss–Hermite rule in each Gaussian variable and an exact sum over the two latent classes.
Endpoints	$d=1$ , $(\theta^{\ast},\nu^{\ast})=(1.10,0.55)$ , and $(\theta,\nu)=(0.45,-0.35)$ . Both the Gaussian covariate and response-noise expectations use 140-point Gauss–Hermite quadrature. For $j=0,\ldots,8$ , the low-noise grid uses $\bar{\alpha}_{t}=1-10^{-6+j/2}$ , while the high-noise grid uses $\bar{\alpha}_{t}=10^{-6+j/2}$ .
Blindness	$d=1$ , $\bar{\alpha}_{t}=\tfrac{1}{2}$ , $(\theta,\nu)=(0.80,-0.35)$ , $\nu^{\ast}=0.55$ , and $\theta^{\ast}=3\times 10^{j/5}>0$ for $j=0,\ldots,25$ . Expectations use deterministic quadrature as detailed below.

Appendix

Table 2: Finite-sample, endpoint, and blindness settings.

Study	Model and numerical setting
Degenerate dynamics	$d=2$ , $c^{\ast}=\mathbf{0}$ , $\eta=1/50=2\times 10^{-2}$ , $\theta_{0}=(\sinh(0.70),0)$ , and $\nu_{0}=0.70$ , so that $J(\vartheta_{0})=0$ . The exact map $\Phi$ is iterated for 30,000 steps and the theorem’s explicit envelope is evaluated at each step.
Nondegenerate dynamics	$d=2$ , $c^{\ast}=(0.80,0)$ , and $\eta=1/50=2\times 10^{-2}$ . Every displayed start is checked to lie in $\mathcal{D}_{\eta}$ before the exact map $\Phi$ is iterated for 1,600 steps. The four starts are $((0.45,0.15),0.40)$ , $((-0.45,0.25),0.40)$ , $((5\times 10^{-3},0.65),0)$ , and the exact saddle $((0,0.65),0)$ .
Spectral diagnostic	$d=20$ , $\Sigma=I_{d}$ , $\theta^{\ast}=(1,0,\ldots,0)$ , $\nu^{\ast}=\nu=0$ , $\bar{\alpha}_{t}=5\times 10^{-3}$ , raw $\theta$ -update step size $\eta_{\theta}=20$ (whose product with $\bar{\alpha}_{t}$ is $10^{-1}$ ), and $\\|\theta_{0}\\|=2\times 10^{-2}$ . The initial absolute cosine with $\theta^{\ast}$ is $10^{-1}$ . For the Gaussian and covariance-normalized spherical Student $t$ designs with ten degrees of freedom, the experiment uses $1.6\times 10^{5}$ common draws and 14 normalized updates.

Appendix

Table 3: Limiting-dynamics and spectral-diagnostic settings.

Sample size $n$	Joint 95% Fisher-ellipse coverage	Relative covariance discrepancy
250	93.5%	0.220
500	96.5%	0.093
1000	94.0%	0.127
2000	94.5%	0.089
4000	95.0%	0.070

Appendix

Table 4: Finite-sample calibration. Joint 95% Fisher-ellipse coverage and relative covariance discrepancy.

Gradient quantity	Low-noise slope (panel a)	High-noise slope (panel b)
$\theta$ component	$1+2.3191\times 10^{-5}$	$1-3.8124\times 10^{-5}$
$\nu$ component	$1+1.0247\times 10^{-5}$	$1-3.00596\times 10^{-4}$

Appendix

Table 5: Endpoint-remainder slopes. Entries are written relative to the unit-slope benchmark.

Figure 5: Auxiliary next-order directional comparison on the balanced slice $\nu^{\ast}=\nu=0$ . Absolute cosine $|\cos\angle(\theta_{k},\theta^{\ast})|$ for normalized finite-noise score-matching updates $\theta\leftarrow\theta-\eta_{\theta}\nabla_{\theta}\mathcal{L}_{t}(\theta,0)$ (solid) and normalized power steps $\theta\mapsto(I_{d}+2\gamma\widehat{\mathcal{C}}^{\ast})\theta/\|(I_{d}+2\gamma\widehat{\mathcal{C}}^{\ast})\theta\|$ with $\gamma=\eta_{\theta}\bar{\alpha}_{t}$ (dotted). Curves labeled $t_{10}$ use spherical Student- $t$ covariates with ten degrees of freedom; the remaining curves use isotropic Gaussians.

Start	$\theta_{0}$	$\nu_{0}$	Role
Same-sign	$(0.45,0.15)$	$0.40$	Attraction to the matching-sign branch.
Opposite-sign	$(-0.45,0.25)$	$0.40$	Crosses toward the target branch.
Near saddle	$(5\times 10^{-3},0.65)$	$0$	Perturbation from the saddle manifold.
Exact saddle	$(0,0.65)$	$0$	Stationary at loss $\\|c^{\ast}\\|^{2}/2$ .

Appendix

Table 6: Nondegenerate-dynamics initializations.

Start	$K_{10^{-2}}$	$K_{10^{-4}}$	$K_{10^{-6}}$	$K_{10^{-8}}$	Tail slope	Tail lack of fit
Same-sign	136	296	475	680	$-0.022705$	$1.71\times 10^{-3}$
Opposite-sign	328	498	694	907	$-0.021873$	$6.67\times 10^{-4}$
Near saddle	614	844	1082	1322	$-0.019248$	$5.68\times 10^{-6}$

Appendix

Table 7: Hitting times and geometric tail fits.

Figure 6: Scope and objective controls for pointwise high-SNR blindness. (a) Fixed-

\bar{\alpha}_{t}

curves provide separate checks covered by the proposition as the true signal

\|\theta^{\ast}\|

grows, for

\bar{\alpha}_{t}\in\{10^{-1},\tfrac{1}{2},1\}

. (b) The moving-candidate curve uses

\theta=0.8/\|\theta^{\ast}\|

, for which the candidate activation need not diverge and the moving-candidate score matching gradient

|\nabla_{\nu}\mathcal{L}_{t}(0.8/\|\theta^{\ast}\|,\nu)|

approaches

0.302302

. That curve is outside the fixed-candidate hypothesis of Proposition 5.2 and illustrates its pointwise-in-candidate scope. (c) The separate population cross-entropy control is evaluated along the same fixed-candidate ray as Figure 4 . It approaches

|\tanh\nu-\tanh\nu^{\ast}|\approx 0.836896

, showing that this selected cross-entropy diagnostic retains imbalance signal.

Score Matching

Score matching estimates gradients of log probability densities to fit models and power diffusion-based generation without computing normalizing constants.

See all 32 papers

Diffusion-based Denoising Beats Vanilla Score Matching in Parameter Estimation: A Theoretical Explanation
May 21, 2026Benedikt Lütke Schwienhorst, Nadja Klein, Johannes LedererScore Matching Parameter Estimation
Department of Mathematics University of Hamburg Hamburg, Germany · Scientific Computing Center Karlsruhe Institute of Technology2026 Karlsruhe, Germany · Department of MathematicsMay University of Hamburg Hamburg, Germany21
Expectation Maximization (EM) Converges for General Agnostic Mixtures
Apr 7, 2026Avishek GhoshExpectation-Maximization Data Mixture Optimization
Department of Computer Science and Engineering, IIT Bombay, Mumbai, India
Generalised Score Matching on Convex Domains
Sep 10, 2026Nishanth Shetty, Saisuchith Mahajan, Chandra Sekhar SeelamantulaScore Matching Parameter Estimation
Department of Electrical Engineering Indian Institute of Science, Bengaluru 560012

May 21, 2026J/K move · Enter open · S save

Diffusion-based Denoising Beats Vanilla Score Matching in Parameter Estimation: A Theoretical Explanation

Benedikt Lütke Schwienhorst, Nadja Klein, Johannes Lederer

Department of Mathematics University of Hamburg Hamburg, Germany · Scientific Computing Center Karlsruhe Institute of Technology2026 Karlsruhe, Germany · Department of MathematicsMay University of Hamburg Hamburg, Germany21

Score Matching Parameter Estimation Score-Based Generative Modeling Diffusion Models

Open paper AskPDF arXiv

Score matching is an alternative to maximum likelihood estimation when the normalizing constant is unknown or too costly to evaluate. However, vanilla score matching has shown to be inefficient relative to maximum likelihood estimation for multimodal distributions with well-separated modes, which are commonly encountered in practical applications. We compare a novel diffusion-based denoising score matching estimator (DDSME) to the vanilla score matching estimator (SME) in this scenario. In particular, we prove statistical guarantees for both estimators, showing that the error bound for the vanilla SME worsens when the separation between the modes increases, which can be avoided in case of the DDSME with suitable hyperparameter tuning. This provides a novel theoretical explanation for the superior behavior of diffusion-based score matching over the vanilla version.

Abstract

Score Matching

Explore similar work

Diffusion-based Denoising Beats Vanilla Score Matching in Parameter Estimation: A Theoretical Explanation

Expectation Maximization (EM) Converges for General Agnostic Mixtures

Generalised Score Matching on Convex Domains