Why Vectorize MAMS Simulation?
A MAMS simulation may require thousands or millions of trial replicates. A simple implementation processes one replicate at a time: generate its data, apply the interim rules, continue or stop, and record the outcome. That approach is easy to inspect, but repeated interpreter-level work can be slow. Vectorization performs the same calculations on arrays of replicate-level inputs, using compiled array operations where available.
The essential constraint is that vectorization must change how calculations are carried out, not what the simulated design does. Each replicate still follows its own path: an arm can be dropped in one replicate and retained in another, and one trial can stop at an interim analysis while another continues. The safe pattern is usually to vectorize across replicates while retaining an explicit loop over the small, fixed number of stages.
Let \(r=1,\ldots,R\) index replicates, \(j=1,\ldots,J\) index stages, and \(a\) index treatment arms. An array might store a stagewise statistic \(Z_{rja}\), with one row per replicate and one column per arm at each stage. A Boolean mask \(M_{rja}\) can indicate whether arm \(a\) is active in replicate \(r\) at stage \(j\). The indices are not merely organizational: keeping replicate, stage, and arm dimensions distinct helps prevent accidental mixing of trial paths.
The data-generating model and its dependence assumptions must be unchanged by the software optimization. For example, if treatment contrasts share a control group, generate the control observations once per replicate and use that same control data in each contrast. Generating an independent control sample for every arm would change the joint distribution of the statistics, even if each individual contrast had the same marginal distribution.
Vectorize Within the Trial’s Stage Structure
Suppose each stage produces a statistic for each arm and applies a prespecified interim rule. The decision at stage \(j\) is a function of the current data and the previous decisions in that replicate. A robust implementation therefore keeps the stage order explicit: calculate stage \(j\) statistics for all currently eligible replicates and arms, update the active masks, and then proceed to stage \(j+1\). A fixed loop over stages is not a failure to vectorize; it preserves the sequential logic while allowing many replicates to be handled together at each stage.
An arm that has been dropped must not become eligible again merely because a later-stage value exists in a pre-generated array. Similarly, a trial that has already stopped must not contribute a later-stage rejection. Masking is the mechanism that carries these rules through vectorized calculations.
Store outcomes or sufficient statistics with replicate and arm indices kept explicit. Preserve shared observations and other required dependence across arms.
Apply the same statistic formula to every eligible replicate, rather than writing a separate calculation for each replicate.
Apply the interim rule separately to each replicate. Carry forward only the arms and trials eligible to continue.
After all applicable decisions are complete, aggregate indicators or other recorded outcomes across replicates.
Worked Example: Vectorizing Stage-One Statistics
Consider three hypothetical replicate snapshots from a two-arm trial with a shared control. At stage one, each arm and control has \(n=25\) observations. Assume independent normal observations within groups, common known standard deviation \(\sigma=10\), and fixed allocation. The control sample mean is generated once per replicate and reused in both treatment contrasts. Suppose the generated means, in outcome units, are:
| Replicate | Control mean | Arm A mean | Arm B mean |
|---|---|---|---|
| 1 | 50 | 54 | 52 |
| 2 | 52 | 53 | 57 |
| 3 | 48 | 51 | 49 |
For either treatment-control contrast, the standard error is the same:
The statistic for arm \(a\) in replicate \(r\) is \(Z_{ra}=(\bar{x}_{ra}-\bar{x}_{rC})/SE\). Applying that formula across all rows gives:
This is a row-wise calculation on arrays: the control-mean column is broadcast across the two arm columns, and the same standard error is used for each contrast. The shared control is not regenerated separately for A and B. If control outcomes are random and the contrasts are correlated because of that shared group, reusing the control values preserves that dependence.
Preserve Dropping and Arm Selection With Masks
Vectorized arithmetic alone does not enforce the design. The decision logic must be represented explicitly. As an illustration, suppose an interim rule retains an arm when its stage-one statistic exceeds \(0\). Among retained arms, select the arm with the largest statistic for a later terminal evaluation. If no arm is retained, the replicate stops without a terminal evaluation. This is only an illustrative rule; an implementation should encode the actual design’s specified thresholds and selection procedure.
For each replicate, form the active-arm mask \(M_{ra}=I(Z_{ra}>0)\). To select the largest statistic without allowing a dropped arm to win, replace inactive values by \(-\infty\) before taking the maximum. The case in which no arm remains must be handled separately: an argmax of two \(-\infty\) values does not mean an arm was selected.
Worked Example: Applying Interim Masks and Selecting an Arm
Suppose four hypothetical replicates have the following stage-one statistics. For each replicate, the retained arms are those with \(Z>0\). Terminal statistics have been generated for both arms, but only the selected arm’s value can determine the terminal decision; use the illustrative terminal boundary \(2.0\).
| Replicate | Stage-one \(Z_A\) | Stage-one \(Z_B\) | Terminal \(Z_A\) | Terminal \(Z_B\) |
|---|---|---|---|---|
| 1 | 1.2 | 0.7 | 2.1 | 1.4 |
| 2 | -0.4 | 1.1 | 2.5 | 1.8 |
| 3 | 0.3 | 0.8 | 1.6 | 2.3 |
| 4 | -0.2 | -0.5 | 2.4 | 2.6 |
In replicate 1, both arms are retained and A is selected because \(1.2>0.7\). Its terminal statistic is \(2.1>2.0\), so the selected arm crosses the illustrative boundary. In replicate 2, only B is retained, so B is selected; \(1.8\) does not exceed \(2.0\). In replicate 3, both arms are retained and B is selected because \(0.8>0.3\); its terminal statistic \(2.3\) exceeds the boundary. In replicate 4, neither arm is retained, so the trial stops and there is no selected-arm terminal decision. The terminal values in that row are ignored.
A vectorized implementation first creates the four-by-two mask. It then assigns inactive stage-one values \(-\infty\), finds the maximum among eligible arms, and evaluates the terminal value only for rows with at least one active arm. The resulting selected-arm decisions are therefore “A, B, B, none,” and the selected arms cross the terminal boundary in replicates 1 and 3 only.
Why the Vectorized and Scalar Calculations Agree
The justification is a coordinate-wise one. Let \(x_r\) contain all realized inputs for replicate \(r\), and let \(F(x_r)\) be the complete deterministic calculation of its statistics and decisions under the specified rules. A scalar implementation evaluates \(F(x_1),\ldots,F(x_R)\) in succession. A correctly written vectorized implementation evaluates those same functions on the corresponding rows. Thus, when both implementations use the same realized inputs and the same row-wise rules, the result for each replicate is the same, except possibly for negligible floating-point differences.
The stronger claim of identical replicate-level results does not follow merely because two data generators have the same joint distribution. If one implementation uses \(U\) and another uses \(1-U\), both inputs can have the same uniform distribution, but the realized values generally differ. If two input arrays have the same joint distribution and the same deterministic decision function is applied, the outputs have the same distribution; they need not match replicate by replicate. This distinction matters when comparing implementations or designs using reproducible random-number streams.
Worked Example: Batch Processing and Replicate Outcomes
Suppose a vectorized MAMS engine processes six hypothetical replicates in batches of two. After applying the same stagewise rules to each row, the indicator of the target event is \(I_r=1\) if that event occurred and \(0\) otherwise. Assume the outputs in replicate order are \(1,0,0,1,0,0\). The first batch contributes \(1+0=1\), the second contributes \(0+1=1\), and the third contributes \(0+0=0\). The total is \(2\) events among \(6\) replicates, giving:
Batching limits memory use because the program need not retain all simulated subject-level outcomes at once. It is valid to add batch event counts and divide by the total replicate count, provided every replicate is included once and the same decision rules are used. If exact replicate-by-replicate comparison with an unbatched run is desired, the batches must use the same realized replicate inputs. Merely using the same input distribution does not guarantee the same six indicators.
Efficiency Choices, Diagnostics, and Limitations
There are three common implementation strategies. A scalar loop over replicates is often easiest to read and debug, but may incur substantial overhead. Full vectorization across replicates can be faster, but it can require large arrays and make branch-specific logic harder to audit. A hybrid approach—vectorized across replicates within each stage, with a loop over stages—is often a practical compromise for MAMS designs because stage order is essential and the number of stages is usually modest.
Vectorization is not automatically faster. Creating huge arrays of individual outcomes for every replicate, arm, and stage can exhaust memory or spend more time moving data than calculating. Chunking the replicate dimension reduces peak memory. Generating later-stage data only for eligible replicates can also reduce unnecessary work, though it requires careful indexing. Pre-generating later-stage values and masking them is simpler in some settings, but inactive values must never influence decisions or summaries.
Useful diagnostics include checking that each replicate has at most one selected arm when the design requires one; that a stopped replicate has no later-stage contribution; and that dropped arms never re-enter. Compare scalar and vectorized implementations on a small set of fixed, realized inputs, including cases with no active arms, one active arm, tied selection statistics, and values near decision boundaries. Near-boundary cases deserve attention because small floating-point differences can change a threshold decision. Define tie handling and numerical tolerances explicitly rather than relying on accidental array ordering.
A common mistake is to apply a terminal decision to every arm in a pre-generated array and then count all boundary crossings, even though the design selected only one arm or stopped earlier. That changes the simulated procedure. The correction is to carry forward replicate-specific eligibility and selection masks, and to count only decisions permitted by the trial’s path. Another mistake is to treat equivalent distributions as identical random-number realizations; distributional agreement is useful, but it is not a pathwise check of the implementation.
Check Your Understanding
Answer each question, then compare with the brief solutions below.
- Why is a loop over stages compatible with vectorization across replicates?
- Two implementations use different random draws but generate inputs with the same joint distribution. What can be concluded about their outputs if they apply the same deterministic decision function?
- In a shared-control design, why should the same simulated control observations be used in both treatment-control contrasts?
- What should a vectorized selection routine do when every arm in a replicate is inactive?
Solutions. (1) Stage order is sequential, while calculations within a stage can still be applied to many replicate rows at once. (2) The output distributions agree under the same deterministic function, but replicate-by-replicate outputs need not be identical. (3) Reusing the control observations preserves the dependence between contrasts induced by their common control group; independent control samples would change the joint distribution. (4) Mark the replicate as having no selected arm and stop its path; do not interpret an argmax among masked values as a valid selection.