Why a Numerical M-Step Is Needed
In SAEM, the E-step’s conditional expectation is replaced by a stochastic approximation built from latent-effect simulation. Some resulting M-steps reduce to formulas, as in a closed-form update. But when a population parameter enters a nonlinear mean function, maximizing the approximated complete-data log likelihood may require numerical optimization. The optimizer is only one part of the update: its target must be the correct SAEM criterion, and its solution must satisfy the model’s parameter constraints.
Consider observations \(y_{ij}\) at times \(t_{ij}\), with subject-specific latent effects \(b_i\) and fixed effects \(\phi\). Write the complete-data log likelihood as \(\ell_c(\phi,\lambda;y,b)\), where \(\lambda\) denotes other model parameters. At SAEM iteration \(k\), the ideal fixed-effect objective, holding \(\lambda\) at its current or separately updated value \(\lambda_k\), is
Here \(\theta_k\) is the current full parameter value and \(\Phi\) is the allowed parameter space for \(\phi\). In SAEM, \(Q_k\) is represented through the current stochastic approximation \(S_k\) of the required complete-data statistics. In an exponential-family representation, this often takes the form \(Q_k(\phi)=-\psi(\phi)+\langle S_k,\nu(\phi)\rangle\), up to terms that do not depend on \(\phi\). When the nonlinear model does not yield a convenient closed-form maximizer, the numerical M-step maximizes this fixed, iteration-specific criterion.
Construct the Criterion Before Calling an Optimizer
For example, suppose a nonlinear population mean is \(f(t;\phi,b_i)\) and residuals are independent normal variables with known variance \(\sigma^2\). The parameter-dependent part of the complete-data log likelihood is
Taking the conditional expectation over \(b\) gives the ideal fixed-effect criterion. In practice, the SAEM statistic or a fixed set of conditional latent-effect draws approximates that expectation. The optimizer must evaluate the same approximation at every candidate \(\phi\). If fresh random draws are generated at each function evaluation, apparent changes in the objective can be simulation noise rather than improvement. Freeze the SAEM approximation for the M-step, or use a numerical scheme explicitly designed for a noisy objective.
Only terms constant in \(\phi\) may be dropped without changing the maximizer. A normalizing term that depends on \(\phi\), for instance, must remain in the criterion. Likewise, maximize \(Q_k\), or equivalently minimize \(-Q_k\); accidentally minimizing \(Q_k\) reverses the intended update. These details are particularly important when the optimizer reports convergence successfully: convergence only describes its numerical progress on the function it was given.
Worked Example: A Scalar Nonlinear Fixed Effect
Suppose the fixed-effect mean is \(f(t;\phi)=e^{-\phi t}\), with \(\phi>0\). For illustration, the current SAEM approximation yields two adjusted responses \(z_1=0.60\) and \(z_2=0.35\), at \(t_1=1\) and \(t_2=2\). These are hypothetical quantities summarizing the current conditional latent-effect calculation; any omitted conditional-variance contribution is constant in \(\phi\). Let \(\sigma^2=0.25\). The criterion to maximize is
Differentiating the displayed criterion gives
A one-dimensional root finder for \(Q'(\phi)=0\), with a bracket over positive values, gives approximately \(\hat{\phi}=0.51905\). At this value, \(e^{-0.51905}\approx0.5950867\) and \(e^{-1.03810}\approx0.3541282\), so
For comparison, at \(\phi=0.50\), the predictions are approximately \(0.6065\) and \(0.3679\), giving \(Q(0.50)\approx-0.000725\). Thus the computed candidate has a larger criterion value. The root is an interior candidate; checking the objective over the allowed domain, including its boundary behavior, is still needed before claiming it is the global maximizer.
Optimization Safeguards and Verification
A numerical M-step should be treated as a constrained maximization problem, not as a request for any stationary point. Positivity constraints can be enforced by optimizing an unconstrained parameter \(u\) with \(\phi=e^u\). This preserves \(\phi>0\), though it can affect numerical scaling. For parameters naturally restricted to an interval, a suitable transformation or a bound-aware optimizer can enforce the interval. In frequentist likelihood optimization, transforming the parameter coordinate does not add a Jacobian term to the likelihood; a Jacobian is relevant when transforming a probability density, not merely when re-expressing its parameter for optimization.
A practical M-step uses sensible starting values, parameter scaling, and a line search or trust region that rejects steps with invalid parameters or a lower objective. If the model has multiple plausible local optima, use several starts and compare the attained values of the same criterion. Then verify the result independently: recalculate the objective, check constraints, inspect the gradient, and assess curvature when feasible. A small gradient indicates approximate stationarity, not global optimality.
Record which SAEM statistics, latent-effect approximation, and nuisance-parameter values are fixed for this M-step. Keep all terms that depend on the target parameter.
Use valid coordinates or explicit bounds, and scale parameters so that differences in units do not make the numerical search unnecessarily ill-conditioned.
Use an optimizer appropriate to the objective; reject invalid or criterion-decreasing trial steps when the method calls for an ascent step. Try additional initial values if local optima are plausible.
Recompute the criterion and gradient at the returned parameter, check constraints and boundary candidates, and compare attained criterion values across starts.
Worked Example: Checking the Direction of Optimization
Use the scalar objective above to check an analytic gradient with a central finite difference at \(\phi=0.50\), using step \(h=0.01\). The finite-difference estimate is
At \(0.51\), the predictions are approximately \(0.6005\) and \(0.3606\), so \(Q(0.51)\approx-[(0.60-0.6005)^2+(0.35-0.3606)^2]/0.50\approx-0.000225\). At \(0.49\), the predictions are approximately \(0.6126\) and \(0.3753\), giving \(Q(0.49)\approx-[(0.60-0.6126)^2+(0.35-0.3753)^2]/0.50\approx-0.00160\). Therefore the finite-difference derivative is
The positive value says that a small increase from \(0.50\) raises \(Q\), consistent with moving toward the maximizer. The analytic derivative gives approximately \(0.06846\) at \(0.50\); the close agreement is a useful implementation check. A sign error that instead treats \(Q\) as a minimization target would move in the wrong direction.
Worked Example: A Two-Parameter Nonlinear Criterion
Consider a hypothetical set of adjusted responses at times \(t=0,1,2\), with values \(2.0000,\ 1.2131,\ 0.7358\). Suppose the nonlinear fixed-effect mean is \(f(t;A,k)=Ae^{-kt}\), with \(A>0\) and \(k>0\), and the residual variance is one. Ignoring constants, the criterion is
At \(A=2,\ k=0.5\), the predictions are \(2.0000,\ 2e^{-0.5}\approx1.2131,\ 2e^{-1}\approx0.7358\). To the displayed precision, the residuals are zero, so \(Q(2,0.5)\approx0\). Since a negative sum of squares cannot exceed zero, this is the global maximum for the unrounded values \(z=(2,2e^{-0.5},2e^{-1})\). With the responses rounded as shown, the residuals are only approximately zero and the maximizer may differ slightly. A numerical check should therefore evaluate the criterion using the actual stored data, not assume that rounded display values give an exact optimum.
For example, the starting point \(A=1.8,\ k=0.4\) predicts \(1.8000,\ 1.2066,\ 0.8088\). Using the displayed responses, the residuals are approximately \(0.2000,\ 0.0065,\ -0.0730\), and \(Q\approx-0.0227\). The fitted candidate has a higher criterion. An optimizer can work with \(u=\log A\) and \(v=\log k\), setting \(A=e^u,\ k=e^v\), so all trial values remain positive. The maximization target remains \(Q(e^u,e^v)\); no Jacobian is added.
Diagnosing Failure and Choosing an Approach
Useful warning signs include a large gradient at the returned point, a criterion value that decreases during an intended ascent step, invalid parameter values, inconsistent results across repeated evaluations of a supposedly frozen objective, or materially different solutions from different starting values. A discrepancy between analytic and finite-difference gradients can indicate a coding or sign error; it can also arise from a poorly chosen difference step or numerical precision, so repeat the check at more than one reasonable step size. A boundary solution requires explicit attention: it may be the constrained maximum, or it may signal that the model or parameter bounds need review.
At an interior local maximum of a twice-differentiable criterion, the Hessian must be negative semidefinite. A negative definite Hessian at a stationary point is sufficient to guarantee a strict local maximum, but it is not necessary; a degenerate strict maximum can have zero curvature in some direction, so higher-order behavior may be needed. Neither a gradient near zero nor favorable local curvature proves global maximality for a non-concave objective.
| Approach | When it is useful | Main caution |
|---|---|---|
| Closed-form M-step | When the current SAEM statistics yield an explicit maximizer. | Do not force a formula when the nonlinear criterion has no closed-form solution. |
| Gradient-based optimization | When reliable derivatives are available and the objective is sufficiently smooth. | May converge to a local optimum; check derivatives, constraints, and starting values. |
| Derivative-free optimization | When derivatives are unavailable or unreliable in a low-dimensional problem. | Can require many evaluations and does not remove the need to verify the objective. |
| Direct Monte Carlo evaluation | When conditional expectations cannot be represented conveniently by stored statistics. | Changing random draws across candidate evaluations can make the target noisy and obscure comparisons. |
A numerical M-step is limited by both the optimizer and the SAEM approximation it receives. A highly accurate maximizer of a noisy or incorrectly constructed criterion is not a reliable parameter update. Conversely, small numerical error may be acceptable when the criterion itself is only approximated, provided the approximation is stable enough for the intended update. The essential distinction is to keep approximation error, optimization error, and model constraints visible rather than treating a convergence flag as proof of correctness.
- Common mistake: trusting the success message. An optimizer can report convergence after minimizing the wrong sign, omitting a parameter-dependent term, or stopping at a poor local solution. Recompute the intended criterion and gradient independently.
- Common mistake: demanding negative definite curvature in every case. Negative definiteness at a stationary point is sufficient for a strict local maximum, not necessary. At an interior local maximum the Hessian is negative semidefinite; degenerate cases may require higher-order checks.
- Common mistake: comparing evaluations with fresh latent draws. Variation may reflect Monte Carlo noise rather than parameter improvement. Hold the SAEM approximation fixed during optimization or account explicitly for a noisy objective.
Check Your Understanding
Answer each question, then compare with the solutions.
- Why should the conditional latent-effect approximation be held fixed while evaluating candidate fixed effects in one numerical M-step?
- If an optimizer minimizes \(Q_k(\phi)\) when the M-step requires maximizing it, what has gone wrong?
- For a positive fixed effect \(\phi\), give a transformation that enforces positivity and state whether a likelihood Jacobian is added.
- At an interior stationary point, what curvature condition is necessary for a local maximum? Is negative definiteness necessary for a strict local maximum?
- A routine reports convergence, but the gradient remains large and results differ across starting values. What should be checked next?
Solutions. (1) Otherwise objective changes can reflect new simulation noise rather than a change in the candidate parameter, making comparisons and the optimization target inconsistent. (2) The direction is reversed; minimizing \(Q_k\) seeks low rather than high values. Minimize \(-Q_k\) instead, or configure the routine to maximize \(Q_k\). (3) Set \(\phi=e^u\) and optimize over real \(u\). No Jacobian is added when this is only a reparameterization of a likelihood optimization. (4) The Hessian must be negative semidefinite. Negative definiteness is sufficient for a strict local maximum at a stationary point, but is not necessary; higher-order behavior may decide a degenerate case. (5) Recompute the objective and gradient, verify the sign and all parameter-dependent terms, check constraints and numerical scaling, and compare valid candidates and attained criterion values across starts.