EM as sufficient-statistic updates

EM separates posterior inference from maximization. A further separation is useful when observations arrive in batches or when estimation is regularized: aggregate posterior statistics, transform that aggregate, then recover model parameters. This mathematical extract of the upstream design develops general and affine update rules and their composition with Bregman shrinkage.

The prerequisites are batch EM, exponential-family coordinates, and mixture architecture. The full upstream framework owns the software interfaces, fitting loops, and product options.

From a batch to a model¶

For a complete-data exponential family with statistic $t(X,Y)$, observation batch $\mathcal B_t$, and current classical parameter $u_{t-1}$, compute

The aggregate contains only the statistic blocks needed by the M-step. The GH convention orders them as $t=(\log Y,Y^{-1},Y,X,X/Y,XX^\top/Y)$. $Y$ is the mixing variable $W$ in the conditioning note; $Z$ remains standard-normal noise. Factor analysis adds latent-factor cross moments to these six blocks.

After an update rule produces a statistic vector $s$, recover parameters by

Here $\theta(u)$ is the natural-coordinate embedding. In an unconstrained minimal family with an attained interior solution, $\theta(\mathcal M(s))=\nabla\phi(s)$. In a curved family, $\mathcal M$ solves the constrained problem; $s$ is not generally the expectation vector of the recovered model. We use the name $\eta$ for the running statistic in the update rules, without asserting ambient moment matching in that case.

Two layers of update rules¶

The general layer is a map with optional state $q_t$:

State can record a cumulative observation count or a running schedule. It is distinct from the sufficient-statistic vector and from the model parameters. A specialization is an affine rule

where $a_t$ is a statistic-space offset and $B_t,C_t$ are linear operators. Scalar weights multiply every block; block-diagonal operators can weight different statistics separately.

In the last row, $n$ is the weight already accumulated and $m$ is the new batch size; the count state becomes $n+m$. The reciprocal schedule $\rho_t=1/(\tau_0+t)$ recovers single-observation weighting with initial weight $\tau_0$. More general decreasing schedules and their assumptions are discussed in online EM.

A scalar convex combination of feasible statistics preserves their convex moment constraints. Arbitrary offsets, operators, or independent block weights need not preserve positive-definiteness, moment inequalities, or existence of an interior maximizing map. The general update notation alone provides no likelihood-ascent or convergence theorem.

Scalar shrinkage as a Bregman penalty¶

Let $\eta_0=\nabla\psi(\theta_0)$ be a reference expectation vector and $\tau\geq0$. The penalized EM derivation maximizes

Dividing by $1+\tau$ makes it ordinary parameter recovery at $S_\tau(\widehat\eta_t)$, where

With a fresh full-data E-step and exact or improving penalized maximization, this is EM for the observed likelihood minus the fixed complete-data KL penalty. In IG dual coordinates, the penalty is $D_\phi(\eta_0\Vert\eta)$; it is not generally a squared Euclidean distance in expectation coordinates.

Composing shrinkage with a running rule¶

Apply the base rule first, then shrink its statistic output:

For an affine base, the composed coefficients are explicit:

Thus one composition handles ordinary EM, forgetting, and sample-weighted updates. If the base is nonlinear, $S_\tau$ is still affine in its own argument, but the complete composed rule need not be affine in its inputs. The base rule’s auxiliary state passes through this construction.

The interpretation depends on the base. With ordinary batch EM, the penalty proof applies to the current observed-data objective. With running statistics, the same algebra describes a regularized surrogate, but does not establish ascent of a fixed observed-data likelihood.

Composition order also matters. Shrinking a running average gives

Updating toward a shrunk new aggregate instead gives

The reference weights are respectively $\tau/(1+\tau)$ and $\rho_t\tau/(1+\tau)$. Repeated shrinkage of the whole running statistic can dominate a vanishing data step: for fixed $\tau>0$, bounded new statistics, and $\rho_t\to0$, the first recurrence tends to $\eta_0$. These schemes therefore need separate choices of schedule and objective.

Blockwise rules and algorithmic limits¶

For nonnegative block weights $\tau_j$, one can write

This defines a blockwise statistic transformation. It is not generally equivalent to the scalar joint-KL penalty, because $\psi$ couples statistic blocks. A separable potential or another explicit penalty derivation would be needed for such a claim. Moment feasibility must also be checked. The covariance-only example explains when changing the sixth GH block alone gives a direct dispersion shrinkage formula.

Batch EM refreshes posterior statistics for all observations at each step. Incremental or mini-batch schemes refresh only the selected observations, retain a running aggregate, and apply $\mathcal M$ to that aggregate. Their fixed points, numerical cost, and convergence depend on the update schedule and on whether old posterior contributions are revisited. Factor constraints, approximate solvers, and post-update transformations must each be assessed against the actual optimization problem.

For numerical inversion and implementation details, see the upstream Bessel and solver design, full EM framework, and package API.

Source and adaptation¶

Rewritten from xshi19/normix, docs/design/em_framework.md, at revision 763bb3608920661a012cf089888d349fbf680aad (2026-09-13 import). Copyright (c) 2020 xshi19. Licensed under MIT. The pinned source records the original version. Notation, mathematical qualifications, and links were adapted for this site. Package interfaces, fitter recipes, and executable cells are omitted; no upstream benchmark or formal-proof verification is claimed.