Maximum Likelihood Iteration

• Maximum likelihood iteration is a family of iterative procedures that maximize or approximate likelihoods in models with latent variables and structural constraints.
• The approach uses alternating updates, gradient-based methods, and fixed-point equations to ensure convergence and feasibility in high-dimensional settings.
• It incorporates approximation techniques and specialized updates to handle singularities and intractable likelihoods in diverse applications such as quantum tomography and state-space models.

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 (Rozzio et al., 9 Mar 2026, Tsai et al., 2022, Lauritzen et al., 2017, Griebel et al., 2019). 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 (Griebel et al., 2019).

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 (Rozzio et al., 9 Mar 2026).

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 (Tsai et al., 2022). A closely related formulation uses outcome matrices $M_n$ and the same negative log-likelihood over $$\theta_\star \in \arg\max_{\theta\in\Theta}\log p_\theta(y), \qquad p_\theta(y)=\int_{\mathcal X} p_\theta(x,y)\,dx.$$0 (Lin et al., 2021). In Gaussian models under total positivity, the MLE takes the log-determinant form

$$\theta_\star \in \arg\max_{\theta\in\Theta}\log p_\theta(y), \qquad p_\theta(y)=\int_{\mathcal X} p_\theta(x,y)\,dx.$$1

maximized over precision matrices $$\theta_\star \in \arg\max_{\theta\in\Theta}\log p_\theta(y), \qquad p_\theta(y)=\int_{\mathcal X} p_\theta(x,y)\,dx.$$2 in the cone of symmetric M-matrices, i.e. $$\theta_\star \in \arg\max_{\theta\in\Theta}\log p_\theta(y), \qquad p_\theta(y)=\int_{\mathcal X} p_\theta(x,y)\,dx.$$3 and $$\theta_\star \in \arg\max_{\theta\in\Theta}\log p_\theta(y), \qquad p_\theta(y)=\int_{\mathcal X} p_\theta(x,y)\,dx.$$4 for $$\theta_\star \in \arg\max_{\theta\in\Theta}\log p_\theta(y), \qquad p_\theta(y)=\int_{\mathcal X} p_\theta(x,y)\,dx.$$5 (Lauritzen et al., 2017).

Other application-specific formulations preserve the same basic pattern. In optical photon-counting, if $$\theta_\star \in \arg\max_{\theta\in\Theta}\log p_\theta(y), \qquad p_\theta(y)=\int_{\mathcal X} p_\theta(x,y)\,dx.$$6 is the probability that a photon lands on pixel $$\theta_\star \in \arg\max_{\theta\in\Theta}\log p_\theta(y), \qquad p_\theta(y)=\int_{\mathcal X} p_\theta(x,y)\,dx.$$7, then the count-pattern log-likelihood is

$$\theta_\star \in \arg\max_{\theta\in\Theta}\log p_\theta(y), \qquad p_\theta(y)=\int_{\mathcal X} p_\theta(x,y)\,dx.$$8

so iterative maximum likelihood becomes repeated maximization of a multinomial-type objective as new photons accumulate (Vella, 2018). In models with analytically intractable likelihoods, the objective itself is replaced by an approximation

$$\theta_\star \in \arg\max_{\theta\in\Theta}\log p_\theta(y), \qquad p_\theta(y)=\int_{\mathcal X} p_\theta(x,y)\,dx.$$9

which is then maximized numerically (Griebel et al., 2019).

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 $$\mathcal F(\theta,q)=\int q(x)\log q(x)\,dx-\int q(x)\log p_\theta(x,y)\,dx,$$0:

$$\mathcal F(\theta,q)=\int q(x)\log q(x)\,dx-\int q(x)\log p_\theta(x,y)\,dx,$$1

The E-step optimizes over probability measures while holding $$\mathcal F(\theta,q)=\int q(x)\log q(x)\,dx-\int q(x)\log p_\theta(x,y)\,dx,$$2 fixed; the M-step optimizes over parameters while holding $$\mathcal F(\theta,q)=\int q(x)\log q(x)\,dx-\int q(x)\log p_\theta(x,y)\,dx,$$3 fixed. In this formulation, the alternation occurs between a Euclidean space for $$\mathcal F(\theta,q)=\int q(x)\log q(x)\,dx-\int q(x)\log p_\theta(x,y)\,dx,$$4 and a Wasserstein space of measures for $$\mathcal F(\theta,q)=\int q(x)\log q(x)\,dx-\int q(x)\log p_\theta(x,y)\,dx,$$5 (Rozzio et al., 9 Mar 2026).

This coordinate-descent interpretation supports particle approximations of the latent distribution. In SVGD-EM, the measure update is represented by interacting particles,

$$\mathcal F(\theta,q)=\int q(x)\log q(x)\,dx-\int q(x)\log p_\theta(x,y)\,dx,$$6

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 (Rozzio et al., 9 Mar 2026).

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 $$\mathcal F(\theta,q)=\int q(x)\log q(x)\,dx-\int q(x)\log p_\theta(x,y)\,dx,$$7 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 (Gong et al., 4 May 2026).

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 (Ramadan et al., 2021).

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,

$$\mathcal F(\theta,q)=\int q(x)\log q(x)\,dx-\int q(x)\log p_\theta(x,y)\,dx,$$8

while the particle update applies the Wasserstein-Nesterov idea in measure space. Algorithmically, the method alternates: update $$\mathcal F(\theta,q)=\int q(x)\log q(x)\,dx-\int q(x)\log p_\theta(x,y)\,dx,$$9 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 (Rozzio et al., 9 Mar 2026).

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

$$\min_{\theta\in\Theta,\;q\in\mathcal P(\mathcal X)} \mathcal F(\theta,q).$$0

and the update

$$\min_{\theta\in\Theta,\;q\in\mathcal P(\mathcal X)} \mathcal F(\theta,q).$$1

Because the Burg entropy keeps iterates in the interior of $$\min_{\theta\in\Theta,\;q\in\mathcal P(\mathcal X)} \mathcal F(\theta,q).$$2, 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 $$\min_{\theta\in\Theta,\;q\in\mathcal P(\mathcal X)} \mathcal F(\theta,q).$$3 and is independent of the sample size $$\min_{\theta\in\Theta,\;q\in\mathcal P(\mathcal X)} \mathcal F(\theta,q).$$4. The expected optimization error of the averaged iterate vanishes at rate

$$\min_{\theta\in\Theta,\;q\in\mathcal P(\mathcal X)} \mathcal F(\theta,q).$$5

(Tsai et al., 2022).

A different first-order construction appears in quantum state tomography via a parameter-free correction of the classical $$\min_{\theta\in\Theta,\;q\in\mathcal P(\mathcal X)} \mathcal F(\theta,q).$$6 method. Starting from $$\min_{\theta\in\Theta,\;q\in\mathcal P(\mathcal X)} \mathcal F(\theta,q).$$7, the update is

$$\min_{\theta\in\Theta,\;q\in\mathcal P(\mathcal X)} \mathcal F(\theta,q).$$8

where

$$\min_{\theta\in\Theta,\;q\in\mathcal P(\mathcal X)} \mathcal F(\theta,q).$$9

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

$\theta$0

Its per-iteration complexity is $\theta$1 (Lin et al., 2021).

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

$\theta$2

with $\theta$3 a damping parameter and $\theta$4 a symmetric negative definite penalization matrix. The paper uses the Marquardt choice $\theta$5 and updates $\theta$6 through a gain ratio comparing the realized increase in log-likelihood to the quadratic prediction. As $\theta$7, the method approaches Newton–Raphson; for large $\theta$8, it becomes more conservative. A local convergence theorem proves linear attraction to the MLE under the stated regularity conditions (Giordan et al., 2014).

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 $\theta$9 and the dual variable $q$0, each updating a low-dimensional block while preserving the M-matrix constraints. Each $q$1-update solves a constrained $q$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 $q$3 observations, irrespective of the underlying dimension (Lauritzen et al., 2017).

A fixed-point form arises in the $q$4-exponential family. There the MLE equation can be written in terms of the dual parameter $q$5 as a convex combination of sufficient statistics with parameter-dependent weights:

$q$6

This yields the iteration

$q$7

For $q$8 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 $q$9 (Tian et al., 6 May 2025).

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 $p_\theta(\cdot\mid y)$0, 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 (Jalali et al., 2022).

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 $p_\theta(\cdot\mid y)$1, a random perturbation $p_\theta(\cdot\mid y)$2 defines

$p_\theta(\cdot\mid y)$3

simulated data are generated under both perturbed parameters, and the log-likelihood gradient is estimated from kernel density estimates on summary statistics:

$p_\theta(\cdot\mid y)$4

The update is then

$p_\theta(\cdot\mid y)$5

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 (Bertl et al., 2015).

The theory of maximum approximated likelihood places such procedures in a general asymptotic framework. MAL defines

$p_\theta(\cdot\mid y)$6

Consistency requires uniform convergence of the approximated criterion to the exact one and a monotonically increasing approximation level $p_\theta(\cdot\mid y)$7. Asymptotic normality requires, in particular, that the gradient approximation error vanish faster than $p_\theta(\cdot\mid y)$8. 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 (Griebel et al., 2019).

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 $p_\theta(\cdot\mid y)$9 physical model applications, and maximum-likelihood computations require $$\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),$$0 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 (Chen et al., 2023).

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

$$\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),$$1

the bound 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),$$2

and therefore

$$\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),$$3

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 (Glancy et al., 2012).

Algebraic statistics studies the global critical-point structure behind such iterations. For a closed irreducible subvariety $$\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),$$4, the maximum likelihood degree $$\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),$$5 is the number of points in the smooth locus where the logarithmic $$\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),$$6-form

$$\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),$$7

degenerates for generic data. In the smooth case,

$$\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),$$8

whereas in the singular case the correct universal bound 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),$$9

The ordinary Euler characteristic does not provide a valid upper bound in general; the paper constructs counterexamples where the gap between $\mathcal D$0 and ML degree is arbitrarily large (Budur et al., 2014).

A related singularity-sensitive object is the maximum likelihood data singular locus $\mathcal D$1, 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

$\mathcal D$2

This identifies a precise mechanism by which the nature of likelihood critical points changes for special data (Horobet et al., 2015). 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.