---
title: Maximum Likelihood Iteration
url: https://www.emergentmind.com/topics/maximum-likelihood-iteration
type: topic
---

# Maximum Likelihood Iteration

Maximum likelihood iteration denotes the class of iterative procedures used to compute maximum-likelihood or maximum-marginal-likelihood estimators when closed-form solutions are unavailable, numerically unstable, or defined on constrained spaces. Across the recent literature, the underlying optimization problem appears in several equivalent forms: direct maximization of a log-likelihood, minimization of a convex negative log-likelihood, alternating optimization over parameters and latent distributions, fixed-point self-consistency equations, and maximization of approximated likelihoods when the exact likelihood is intractable [2603.08676] [2211.12880] [1702.04031] [1908.04110]. This suggests that maximum likelihood iteration is best understood not as a single algorithm, but as a family of update principles adapted to latent variables, matrix constraints, non-Euclidean geometry, singular models, and large-scale numerical linear algebra.

## 1. Objective structure and problem classes

In its standard form, maximum likelihood estimation seeks
$$
\hat\theta \in \arg\max_{\theta\in\Theta} \frac{1}{n}\sum_{i=1}^n \log f(z_i;\theta),
$$
or equivalently minimizes the negative log-likelihood. Several recent formulations emphasize that the computational character of the iteration is determined less by the likelihood principle itself than by the geometry of the feasible set and by whether latent variables, implicit integrals, or structural constraints are present [1908.04110].

For latent-variable models, a central example is maximum marginal likelihood estimation (MMLE), where
$$
\theta_\star \in \arg\max_{\theta\in\Theta}\log p_\theta(y), \qquad
p_\theta(y)=\int_{\mathcal X} p_\theta(x,y)\,dx.
$$
The free-energy reformulation introduces
$$
\mathcal F(\theta,q)=\int q(x)\log q(x)\,dx-\int q(x)\log p_\theta(x,y)\,dx,
$$
and rewrites MMLE as a joint minimization over parameters and probability measures,
$$
\min_{\theta\in\Theta,\;q\in\mathcal P(\mathcal X)} \mathcal F(\theta,q).
$$
For fixed $\theta$, minimization over $q$ yields the posterior $p_\theta(\cdot\mid y)$, and the MMLE optimum corresponds to the global optimum of the free-energy problem [2603.08676].

In constrained matrix estimation, the likelihood objective can remain convex while the domain is highly structured. In maximum-likelihood quantum state tomography, the estimator is
$$
\hat{\rho}\in\arg\min_{\rho\in\mathcal D} f(\rho),\qquad
f(\rho)=\frac{1}{n}\sum_{i=1}^n -\log \operatorname{tr}(A_i\rho),
$$
with $\mathcal D$ the set of density matrices [2211.12880]. A closely related formulation uses outcome matrices $M_n$ and the same negative log-likelihood over $\mathcal D$ [2110.00747]. In Gaussian models under total positivity, the MLE takes the log-determinant form
$$
\ell(K;S)=\log\det K-\operatorname{tr}(SK),
$$
maximized over precision matrices $K$ in the cone of symmetric M-matrices, i.e. $K\succ 0$ and $K_{ij}\le 0$ for $i\neq j$ [1702.04031].

Other application-specific formulations preserve the same basic pattern. In optical photon-counting, if $P(i\mid \mathbf p)$ is the probability that a photon lands on pixel $i$, then the count-pattern log-likelihood is
$$
\ell(\mathbf p\mid \mathbf I)=\ln P_0+\sum_i I_i \ln P(i\mid \mathbf p),
$$
so iterative maximum likelihood becomes repeated maximization of a multinomial-type objective as new photons accumulate [1806.04503]. In models with analytically intractable likelihoods, the objective itself is replaced by an approximation
$$
\tilde L_n(\theta)=\frac{1}{n}\sum_{i=1}^n \log \tilde f_{R(n)}(z_i;\theta),
$$
which is then maximized numerically [1908.04110].

## 2. Alternating and latent-variable iterations

Expectation-maximization is the canonical alternating maximum-likelihood iteration. In the free-energy formulation of MMLE, EM becomes coordinate descent on the joint variable $(\theta,q)$:
$$
q_t=\arg\min_q \mathcal F(\theta_t,q),\qquad
\theta_{t+1}=\arg\min_\theta \mathcal F(\theta,q_t).
$$
The E-step optimizes over probability measures while holding $\theta$ fixed; the M-step optimizes over parameters while holding $q$ fixed. In this formulation, the alternation occurs between a Euclidean space for $\theta$ and a Wasserstein space of measures for $q$ [2603.08676].

This coordinate-descent interpretation supports particle approximations of the latent distribution. In SVGD-EM, the measure update is represented by interacting particles,
$$
q_t \approx \frac{1}{N}\sum_{i=1}^N \delta_{x_t^{(i)}},
$$
and evolved deterministically through Stein variational gradient descent rather than Langevin dynamics. The resulting particle evolution is deterministic, interacting or nonlocal through the kernel, and intended as a practical approximation to free-energy descent [2603.08676].

Several EM generalizations refine the alternation by changing the complete-data representation. In deterministic maximum-likelihood direction finding with Gaussian mixture noise, SAGE uses multiple EM-pairs, whereas AECM employs multiple less informative complete-data versions and updates the directions of arrival sequentially, one by one, rather than simultaneously. The AECM construction contains $M+2$ EM-cycles per iteration, uses golden section search to locate the local maximum closest to the previous estimate, and is shown to converge to a stationary point under standard regularity conditions. The numerical results reported there indicate faster stable convergence than SAGE, with almost the same per-iteration complexity [2605.02309].

The same EM logic has been adapted to recursive state estimation. For nonlinear and non-Gaussian state-space models, a particle filter approximates the conditional density, and the EM algorithm is then used to maximize the posterior density with respect to the current state. The filtering, prediction, and smoothing variants all inherit the standard EM ascent property: the posterior likelihood sequence is monotonically non-decreasing, and under unimodality the iteration converges to the unique stationary point. With randomized reinitialization, the method is reported to reach the MLE for multimodal, truncated, skewed, and disjoint-support densities [2103.10475].

## 3. First-order, mirror, and accelerated schemes

A major development in maximum likelihood iteration is the replacement of exact or implicit M-steps by first-order or particle-based updates with explicit complexity guarantees. In MMLE, the Euclidean–Wasserstein gradient-flow viewpoint yields alternating parameter and measure updates, and the paper on Momentum SVGD-EM augments both with Nesterov acceleration. The parameter update uses the classical look-ahead mechanism,
$$
\theta_{t+1}=\tilde\theta_t+\frac{\gamma}{N}\sum_{i=1}^N \nabla_\theta \ell(\tilde\theta_t,x_t^{(i)}),\qquad
\tilde\theta_{t+1}=\theta_{t+1}+\alpha_\theta(\theta_{t+1}-\theta_t),
$$
while the particle update applies the Wasserstein-Nesterov idea in measure space. Algorithmically, the method alternates: update $\theta$ with Nesterov momentum, update particles with SVGD, apply particle momentum, and repeat. The reported experiments show consistent acceleration in iteration count, with particle-space acceleration identified as the main driver of speedup as dimension grows [2603.08676].

In quantum tomography, first-order likelihood iteration must respect positivity and trace constraints while avoiding singular iterates. Stochastic mirror descent with the Burg entropy uses the mirror map
$$
h(\rho)=-\log\det \rho,
$$
and the update
$$
\rho_{t+1}\leftarrow \arg\min_{\rho\in\mathcal D}
\left\{
\eta\, \operatorname{tr}\!\big(\nabla f_t(\overline\rho_t)(\rho-\rho_t)\big)
+ D_h(\rho,\rho_t)
\right\}.
$$
Because the Burg entropy keeps iterates in the interior of $\mathcal D$, it avoids the singularity and stalling issues that can break projected gradient descent. The bottleneck is a single eigendecomposition, so the per-iteration time is $O(d^3)$ and is independent of the sample size $n$. The expected optimization error of the averaged iterate vanishes at rate
$$
O\!\left(\sqrt{\frac{1}{T}d\log T}\right)
$$
[2211.12880].

A different first-order construction appears in quantum state tomography via a parameter-free correction of the classical $R\rho R$ method. Starting from $\rho_1=I/D$, the update is
$$
\rho_{k+1}=\mathcal N\!\left(\exp\!\left(\log(\rho_k)+\log(R(\rho_k))\right)\right),
$$
where
$$
R(\rho_k)=\frac{1}{N}\sum_{n=1}^N \frac{M_n}{\operatorname{tr}(M_n\rho_k)}.
$$
The method is full-rank, needs no step-size tuning, reduces to Cover’s multiplicative update in the commuting case, and satisfies the non-asymptotic bound
$$
f(\overline{\rho}_k)-f(\hat{\rho})\le \frac{\log D}{k}.
$$
Its per-iteration complexity is $O(D^3+ND^2)$ [2110.00747].

## 4. Constraint-preserving, damped, and fixed-point updates

Many maximum-likelihood iterations are designed to preserve feasibility and improve conditioning rather than to minimize iteration count alone. In exponential-family models, a Levenberg–Marquardt adaptation replaces the raw Newton step by
$$
\beta^{(t+1)}=\beta^{(t)}-\big[H_l(\beta^{(t)})+\gamma^{(t)}P(\beta^{(t)})\big]^{-1}s(\beta^{(t)}),
$$
with $\gamma^{(t)}>0$ a damping parameter and $P(\cdot)$ a symmetric negative definite penalization matrix. The paper uses the Marquardt choice $P(\theta^{(t)})=\operatorname{diag} H_l(\theta^{(t)})$ and updates $\gamma$ through a gain ratio comparing the realized increase in log-likelihood to the quadratic prediction. As $\gamma^{(t)}\to 0$, the method approaches Newton–Raphson; for large $\gamma^{(t)}$, it becomes more conservative. A local convergence theorem proves linear attraction to the MLE under the stated regularity conditions [1410.0793].

In Gaussian MTP2 models, feasibility is defined by sign constraints on the precision matrix, and the resulting maximum-likelihood iteration is coordinate descent. The paper gives globally convergent coordinate-descent algorithms on both the primal variable $K$ and the dual variable $\Sigma$, each updating a low-dimensional block while preserving the M-matrix constraints. Each $K$-update solves a constrained $2\times 2$ maximization subproblem; the method is explicitly described as structurally similar to iterative proportional scaling. The same work proves that the MLE exists with probability one for at least $2$ observations, irrespective of the underlying dimension [1702.04031].

A fixed-point form arises in the $\lambda$-exponential family. There the MLE equation can be written in terms of the dual parameter $\eta$ as a convex combination of sufficient statistics with parameter-dependent weights:
$$
\hat\eta=\sum_{i=1}^n w_i(\hat\theta)\,y_i.
$$
This yields the iteration
$$
\eta(k+1)=\sum_{i=1}^n w_i(\theta(k))\,y_i,\qquad
\theta(k+1)=\nabla^{(\lambda)}\psi(\eta(k+1)).
$$
For $\lambda<0$ and under the paper’s regularity assumptions, the likelihood strictly increases along the iterates unless the current point is already a fixed point. The monotonicity proof is based not on EM, but on the strict concavity of $\Psi(\eta)=e^{\lambda\psi(\eta)}$ [2505.03582].

An online adaptive-filtering variant appears in ML-APA, where maximum likelihood is used to adapt the regularization in the affine projection algorithm. Under memoryless Gaussian inputs, the expected misalignment converges to zero as $O(1/t)$, and the achieved error is asymptotically optimal. The incremental version, IML-APA, updates at every time step and outperforms ML-APA in the reported simulations [2209.01594].

## 5. Approximate and scalable likelihood iteration

When exact likelihood evaluation is analytically unavailable or computationally prohibitive, maximum likelihood iteration shifts from exact optimization to optimization of an approximation. In the approximate maximum likelihood method for intractable models, the iteration uses simultaneous perturbation stochastic approximation. At step $k$, a random perturbation $\delta_k$ defines
$$
\Theta^\pm=\Theta_{k-1}\pm c_k\delta_k,
$$
simulated data are generated under both perturbed parameters, and the log-likelihood gradient is estimated from kernel density estimates on summary statistics:
$$
\hat{\nabla}_{c_k} \hat{l}(\Theta_{k-1}; S_{obs})
=
\delta_k \frac{\log \hat{L}(\Theta^+; S_{obs}) - \log \hat{L}(\Theta^-; S_{obs})}{2c_k}.
$$
The update is then
$$
\Theta_k=\Theta_{k-1}+a_k \hat{\nabla}_{c_k}\hat l(\Theta_{k-1};S_{obs}).
$$
The method is local rather than global, differs from ABC by following a simulated gradient instead of sampling broadly across the parameter space, and is explicitly recommended to be run from multiple random starting values to reduce the risk of local maxima [1507.04553].

The theory of maximum approximated likelihood places such procedures in a general asymptotic framework. MAL defines
$$
\hat\theta^{(n)}_{MAL}
=
\arg\max_{\theta\in\Theta}\tilde L_n(\theta),\qquad
\tilde L_n(\theta)=\frac{1}{n}\sum_{i=1}^n \log \tilde f_{R(n)}(z_i;\theta).
$$
Consistency requires uniform convergence of the approximated criterion to the exact one and a monotonically increasing approximation level $R(n)\to\infty$. Asymptotic normality requires, in particular, that the gradient approximation error vanish faster than $n^{-1/2}$. The framework covers Monte Carlo, quasi-Monte Carlo, Gaussian quadrature, and sparse grids, and the paper’s comparison shows that deterministic high-order quadrature can require far fewer points than Monte Carlo in smooth problems [1908.04110].

At larger scales, the iteration can remain exact at the optimization level while using compressed linear algebra inside each likelihood evaluation. In physics-based Gaussian process models, hierarchical off-diagonal low-rank approximations replace dense covariance matrices and their derivatives. The construction requires $O(\log n)$ physical model applications, and maximum-likelihood computations require $O(n\log^2 n)$ effort per iteration. The same HODLR structure is differentiated to approximate covariance derivatives, and exact trace formulas give the score and expected Fisher information in the same quasilinear complexity class [2303.10102].

## 6. Stopping rules, critical-point geometry, and singularities

Maximum likelihood iteration is governed not only by update formulas but also by criteria for termination and by the geometry of the critical points being sought. In quantum tomography, a gradient-based upper bound controls the residual optimality gap independently of the particular iterative algorithm. With
$$
R(\rho_k)=\sum_{i=1}^{N} \frac{\Pi_i}{\Tr(\rho_k\Pi_i)},
$$
the bound is
$$
L(\rho_{\mathrm{ML}})-L(\rho_k)\le r_k,\qquad
r_k=\max\{\mathrm{eig}[R(\rho_k)]\}-N,
$$
and therefore
$$
\frac{\mathcal L(\rho_{\mathrm{ML}})}{\mathcal L(\rho_k)}\le e^{r_k}.
$$
This yields stopping rules for point estimation, likelihood-ratio confidence regions, and expectation-value confidence intervals, with Wilks’s theorem providing the relevant chi-squared scales [1205.4043].

Algebraic statistics studies the global critical-point structure behind such iterations. For a closed irreducible subvariety $X\subset(\mathbb C^*)^n$, the maximum likelihood degree $\mathrm{MLdeg}(X)$ is the number of points in the smooth locus where the logarithmic $1$-form
$$
\lambda_1 \frac{dp_1}{p_1}+\cdots+\lambda_n \frac{dp_n}{p_n}
$$
degenerates for generic data. In the smooth case,
$$
(-1)^d\chi(X)=\mathrm{MLdeg}(X),
$$
whereas in the singular case the correct universal bound is
$$
(-1)^d \chi\bigl(IC(\mathbb C_X)\bigr)\ge \mathrm{MLdeg}(X).
$$
The ordinary Euler characteristic does not provide a valid upper bound in general; the paper constructs counterexamples where the gap between $\chi$ and ML degree is arbitrarily large [1411.3486].

A related singularity-sensitive object is the maximum likelihood data singular locus $\mathrm{DS}(X)$, the set of data for which a likelihood critical point lies in the singular locus of the model away from the coordinate hyperplanes. The paper proves the Hadamard-product bounds
$$
(X_{sing}\setminus H)*[1:\ldots:1:-1] \subseteq \mathrm{DS}(X)
\subseteq (X_{sing}\setminus H)*X^*.
$$
This identifies a precise mechanism by which the nature of likelihood critical points changes for special data [1509.09225]. This suggests that, in singular models, the behavior of maximum likelihood iteration cannot be understood solely through local descent or ascent properties; it also depends on how the likelihood equations interact with singular strata and with the global count of critical points.

Source: https://www.emergentmind.com/topics/maximum-likelihood-iteration