What Can Run at the Same Time?
A SAEM iteration contains both independent computations and computations that depend on a shared, current result. The practical opportunity for parallelism is usually the first kind: conditional simulations and statistic calculations for different subjects. The main constraint is that the stochastic-approximation update must use the complete statistic for the iteration, not a partial sum that happens to be available first.
Let \(y_i\) denote the observations for subject \(i\), \(b_i\) that subject’s latent effects, and \(\theta\) the shared population-parameter vector. There are \(N\) subjects. Write \(t_i(b_i,y_i)\) for subject \(i\)’s contribution to the vector of complete-data sufficient statistics, and let \(T\) be the corresponding total. In the usual subject-level setup, the complete-data log likelihood is a sum over subjects, and the statistic contributions can therefore be added.
This factorization is an assumption to check, not a consequence of using a computer with multiple processors. It is typical when subjects have independent latent effects conditional on shared population parameters and their own observations. If the model includes a shared latent effect linking subjects, then the conditional draws may be dependent. In that case, simulating subjects independently generally samples the wrong joint distribution.
For iteration \(k\), suppose the latent effects are simulated using the current parameter vector \(\theta_{k-1}\). Once all subject contributions for that iteration are available, define \(T_k^{\mathrm{sim}}=\sum_{i=1}^N t_i(b_i^{(k)},y_i)\). A common stochastic-approximation update is
where \(S_{k-1}\) is the stored statistic vector and \(\gamma_k\) is the iteration’s step size. The M-step then uses \(S_k\) to update the shared parameter vector. This tutorial uses a total statistic; an implementation based on average per-subject statistics must use the matching normalization consistently.
Why the Synchronization Barrier Matters
Subject-level simulations and evaluations can be assigned to different workers because each task uses the same iteration-level parameters and its own subject’s data. Each worker returns its statistic contribution, or a partial sum of contributions. Before applying the stochastic-approximation update, the program must wait until every required contribution for that iteration has arrived and must combine them into \(T_k^{\mathrm{sim}}\). That wait is a synchronization barrier.
The barrier is necessary because \(S_k\) is defined using the full sum for iteration \(k\). Updating \(S_k\) as soon as only some subjects finish changes the recursion: the result then depends on task completion order or omits contributions. The parameter update must in turn wait for \(S_k\), since it is computed from that updated statistic. Thus the broad dependency is: simulate and calculate subject contributions in parallel; combine them; update \(S_k\); perform the shared M-step; then begin computations that require the new \(\theta_k\).
The justification is the additivity of the statistic. If the subjects are divided among \(R\) workers, with worker \(r\) assigned the index set \(I_r\), let \(T_{k,r}^{\mathrm{sim}}=\sum_{i\in I_r}t_i(b_i^{(k)},y_i)\). The sets \(I_r\) partition the subjects, so
Substituting this same total into the SA recursion yields the same \(S_k\). If the M-step is also applied to the same \(S_k\), it receives the same input as in the serial calculation. This argument concerns the calculation conditional on the same draws; in practical programs, independent random-number streams can generate different draws, and floating-point summation can produce tiny numerical differences.
Worked Example: Reducing Contributions Before the SA Update
Suppose a hypothetical iteration has four subjects. Their simulated contributions to one sufficient statistic are 4, 8, 6, and 10 units. The previous statistic is \(S_{k-1}=10\) units, and the step size is \(\gamma_k=0.25\). Assume each task uses the same current parameter vector and the contributions are additive.
Two workers can calculate the first two and last two contributions, respectively. Their partial sums are \(4+8=12\) units and \(6+10=16\) units, so the complete simulated statistic is \(T_k^{\mathrm{sim}}=12+16=28\) units.
The synchronized update is
Thus the updated statistic is 14.5 units. The workers can form their partial sums concurrently, but the update is not ready until both sums have been combined. Applying the update to only the first worker’s partial sum would instead give \(0.75(10)+0.25(12)=10.5\) units, which is not the specified iteration update.
Parallel Work and Sequential Dependencies
The most straightforward distribution is across subjects. A worker can simulate the current conditional latent-effect state for an assigned subject, evaluate that subject’s contribution, and return the contribution to the reduction. Subject tasks may differ in duration, so a practical scheduler can assign new tasks to workers as they become available. This can reduce idle time, but it does not change the requirement to collect every contribution before updating the shared statistic.
For a subject whose conditional simulation uses a Markov chain, successive transitions within that chain ordinarily depend on the preceding state. Those transitions are not made independent merely by splitting their computation across workers. Across-subject parallelism remains available when the subject-level conditional distributions factorize, but the chain evolution for a given subject must respect its transition dependence. Parallelizing independent chains is possible in some designs, but those chains must represent the intended simulation procedure; they do not justify treating dependent successive states as independent draws.
The shared M-step is another boundary. If it maximizes an objective defined from the global updated statistic, it must use the fully reduced \(S_k\). For a numerical M-step, evaluations of subject-level objective or derivative contributions may sometimes be computed in parallel, but the global evaluation must combine the contributions before the optimizer takes a step. Whether a particular parameter update separates into independent pieces depends on the model and constraints; shared parameters or covariance constraints can couple components.
Make the current parameter vector and any iteration-specific settings available to all subject tasks. Do not let workers silently use different parameter versions.
Simulate each subject’s latent effects under the intended conditional distribution and calculate that subject’s statistic contribution.
Confirm that every required task has completed, then sum the subject contributions to form \(T_k^{\mathrm{sim}}\).
Apply the SA recursion to the complete statistic, carry out the M-step, and make the resulting parameters available to the next iteration.
Worked Example: Estimating the Benefit of Subject-Level Parallelism
In a hypothetical run with 12 subjects, suppose each subject’s conditional simulation and statistic calculation takes 8 seconds. Ignore variation in task times for this calculation. A serial iteration spends \(12(8)=96\) seconds on those tasks. With four workers, each receiving three subjects, the subject tasks take \(3(8)=24\) seconds if all workers stay busy. Suppose collecting and reducing the results takes an additional 3 seconds.
The parallel iteration therefore takes approximately \(24+3=27\) seconds, and the speedup for this part of the iteration is \(96/27=3.56\). This is less than the ideal fourfold speedup because the reduction adds time. The calculation also assumes the tasks balance evenly; if one worker receives a much slower subject, other workers may wait for it at the barrier. These figures describe a hypothetical computational setup, not a guarantee for a particular model or machine.
If other serial parts of an iteration take substantial time, the total iteration speedup will be smaller still. Parallelizing subject tasks cannot make a serial shared-parameter update disappear.
Checks, Comparisons, and Limitations
Compare parallel and serial calculations on a small test problem using the same simulated draws when possible. Check that every subject contributes exactly once, the reduced total matches a direct sum, the SA update uses the intended step size, and the M-step starts only after that update. Small discrepancies can arise because floating-point addition is not exactly associative: different reduction orders may round \( (a+b)+c \) and \( a+(b+c) \) differently. Large or systematic discrepancies suggest a more consequential issue, such as omitted tasks, duplicate contributions, inconsistent parameter versions, or a reduction performed at the wrong time.
Parallel subject-level computation is most effective when subject tasks are sufficiently costly, reasonably balanced, and conditionally independent. For small tasks, communication and scheduling overhead can exceed the saved computation. With distributed workers, transferring latent states or partial statistics can also be costly. A practical implementation should record which subject tasks completed and should fail clearly if a required contribution is missing rather than updating from an incomplete sum.
| Approach | What it does | Main consideration |
|---|---|---|
| Serial computation | Processes subject tasks one at a time, then updates the global statistic. | Simple to inspect; may be slow when many subject tasks are costly. |
| Synchronized parallel reduction | Processes eligible subject tasks concurrently, then combines all contributions before updating. | Preserves the specified iteration calculation for the same draws, subject to numerical rounding. |
| Asynchronous shared updates | Allows workers to update shared statistics as their tasks finish. | Can use partial or stale values and make the update order-dependent; it is not the same recursion without a separately justified algorithm. |
The synchronous approach is not always the fastest possible design, but it has a clear statistical interpretation: it computes the complete iteration statistic before applying the stated SA update. An asynchronous method may be appropriate only if its modified update rule and convergence properties are justified for the problem; it should not be treated as an implementation-only change. The reproducibility of stochastic SAEM runs also benefits from recording how tasks are assigned and how random-number streams are managed.
Common Mistake
This pitfall is not fixed by averaging the worker-specific updates: averaging 3.5 and 5.5 gives 4.5, also not 6.5. The update is defined using the full total contribution, so the reduction must occur before the single global update.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions below.
- What conditional-independence condition supports parallel simulation of subject-level latent effects?
- For \(S_{k-1}=8\), \(\gamma_k=0.2\), and three contributions of 5, 7, and 9 units, calculate \(S_k\).
- Why can’t an iteration use the new \(\theta_k\) for some subjects while other subjects are still being simulated using \(\theta_{k-1}\), unless a different algorithm has been specified?
- Name one reason that a parallel iteration may be slower than its ideal worker-count speedup suggests.
Solutions. (1) The relevant condition is factorization of the conditional distribution, \(p(b_1,\ldots,b_N\mid y_1,\ldots,y_N,\theta)=\prod_i p(b_i\mid y_i,\theta)\), along with subject-specific tasks that do not require another subject’s draw. (2) The total is \(5+7+9=21\) units, so \(S_k=(1-0.2)(8)+0.2(21)=6.4+4.2=10.6\) units. (3) The iteration’s contributions are defined using one shared parameter vector; mixing parameter versions changes the calculation and no longer gives the stated SAEM iteration. (4) Examples include synchronization and communication overhead, uneven task durations, or a serial M-step.