Computer Science Computer Science major 100% Free Open Access
Chapter 7 • Theory & Derivations

Numerical and data workloads

Simulation; Graph Processing; Machine Learning; Scientific Pipelines

§7.1 Simulation

Numerical simulation transforms equations into discrete updates that a computer can execute. Many physical models consist of local updates, making them promising for parallelization, but the exact dependencies depend on the numerical method. An explicit finite-difference scheme for heat diffusion calculates the next time step from the current grid state. Each cell's update reads neighboring values but need not write to another cell's output location. Separate old and new arrays give a clean parallel ownership rule. An implicit scheme instead requires solving a coupled linear system at each time step, introducing different synchronization and communication challenges. Parallelism follows the mathematical discretization, not the word 'simulation' alone.

Consider a one-dimensional heat equation discretized as u_i^(k+1)=u_i^k+r(u_(i-1)^k-2u_i^k+u_(i+1)^k), where r is the diffusion coefficient times time step divided by squared spatial step. For the standard explicit scheme with simple boundaries, stability requires r not exceed one-half under familiar assumptions. Assign contiguous index ranges to p workers, and let each compute output values from the unchanged old array. At rank boundaries, exchange one ghost value from each neighbor before updating the outermost cells. After all outputs are complete, swap old and new buffers. A worker that updates the old array in place violates the scheme's intended time-level semantics and may cause schedule-dependent results even if it never writes another worker's owned index.

Communication costs arise from subdomain boundaries. For a three-dimensional stencil, each rank may exchange six faces, and corner exchanges may be necessary for more complex operators. The local compute cost scales with volume, while exchanged cells scale more closely with surface area. As the global problem is spread over more ranks, the surface-to-volume ratio increases and strong-scaling efficiency can fall. Increasing stencil order also widens halos. Some algorithms perform several local time steps between exchanges by temporarily extending halos, trading redundant arithmetic and storage for fewer communication startups. The numerical equivalence of such methods must be established rather than assumed.

Monte Carlo simulation offers a different model. Independent trajectories can often run with little communication until their aggregate statistics are reduced. The main correctness concern becomes random-number stream independence and reproducibility. Assigning each worker an identical seed can produce perfectly correlated trajectories, invalidating the estimated uncertainty. Simply changing a seed by one does not constitute a proof of independence for an arbitrary generator. Counter-based generators or carefully partitioned streams can provide well-defined reproducible assignment. If results are aggregated in a different floating-point order, small numerical differences may still occur; record the reduction method.

Adaptive simulations produce irregular work. An adaptive mesh refines some regions far more than others, so assigning equal numbers of coarse cells may not balance computational cost. A fluid solver may spend extra time where turbulence is intense; a particle code may cluster many particles into one spatial region. Load balancing then becomes an ongoing process. Dynamic repartitioning requires transferring physical state, rebuilding connectivity and ensuring no cells or particles are lost. An important decision is whether the saved cost per future iteration exceeds the one-time migration cost. More frequent balancing is not automatically better.

Checkpointing matters because long parallel simulations can run for hours or days and use thousands of components. Saving state periodically protects against failures but consumes I/O bandwidth and storage. If each checkpoint takes considerable time, checkpointing too often reduces useful throughput; too infrequent checkpointing risks losing large amounts of progress. The choice depends on failure rates, checkpoint duration and recovery cost. Consistent state requires that checkpoint files represent a coherent simulation time level. Partial files must not be mistaken for completed checkpoints, and metadata should include dimensions, step count, numerical scheme and version information.

Numerical validation cannot rely only on comparing one large parallel run with another. For heat diffusion, test constant fields, for which the update should remain constant under compatible boundary conditions. Test conservation where the model implies it, and compare against analytical solutions or high-resolution references. Run small problems with one and multiple workers and compare within justified tolerances. Use a convergence study as spatial and temporal resolution change; correct parallelization does not rescue an unstable or inaccurate discretization. A highly scalable simulation that solves the wrong equations is an engineering failure.

Performance reports should name the simulated physical time, numerical resolution, time-step count, solver tolerance, hardware and initialization costs. Weak scaling can answer whether more hardware supports higher resolution within the same wall time, while strong scaling measures time saved for one fixed scientific problem. The final goal is reliable science, not merely a good speedup curve.

§7.2 Graph Processing

Graphs represent irregular relationships through vertices and edges. Parallel graph computation is difficult because degrees and neighborhoods vary, memory accesses are often indirect, and the active work can change dramatically between iterations. Breadth-first search, connected components, PageRank and graph neural network message aggregation all use graph structure but expose different opportunities for parallelism. A simple vertex-parallel implementation assigns one vertex to each worker, yet a vertex with a million edges performs far more work than one with two edges. Balancing vertex counts does not therefore balance execution time.

Breadth-first search illustrates dynamic frontiers. Start with one source vertex at distance zero. At each level, inspect edges from the current frontier and identify previously unvisited neighbors for the next frontier. Every newly discovered vertex must receive its correct distance and must not be processed in conflicting ways. Multiple workers may discover the same neighbor concurrently, so visited-state updates require an atomic claim, a partitioned owner, or another deduplication mechanism. Duplicate discoveries that are not handled correctly can waste work or corrupt output. A level-synchronous algorithm uses barriers between frontiers; an asynchronous alternative requires more elaborate reasoning about ordering and termination.

A traditional compressed sparse row representation stores an offset for each vertex and contiguous arrays of neighboring vertex IDs. This format avoids pointer-heavy adjacency objects and often improves sequential edge traversal. Yet neighbor IDs lead to irregular reads of visited flags and vertex state, which can miss caches. GPU traversal may process many edges per active vertex cooperatively to reduce imbalance. On a CPU, a task scheduler might assign chunks of edges rather than complete vertices. The representation should be chosen by the algorithm's access pattern and the expected graph size, not merely by familiarity with one library's API.

Graph partitioning for a cluster aims to assign vertices or edges to ranks while limiting communication. A partition with very few cut edges can be attractive, but it may place most high-degree vertices on one rank, causing imbalance. Vertex cuts replicate high-degree vertex state across partitions so that edges can be balanced, at the cost of synchronization and increased state replication. Edge cuts keep vertices uniquely owned and communicate across cut edges. Neither strategy dominates every graph. Social networks with heavy-tailed degree distributions behave differently from regular mesh graphs or road networks.

PageRank-type iterative computations propagate values along graph edges and then combine contributions. Their parallel implementation resembles a sparse matrix-vector multiplication or a message-passing update. A 'push' approach sends a vertex's contribution to its neighbors, potentially causing conflicting writes to destination accumulators. A 'pull' approach has each destination gather contributions from incoming neighbors, making output ownership easier but changing memory-access behavior. The best choice depends on graph density, frontier size and data layout. Hybrid algorithms switch modes based on the active fraction of vertices.

Distributed graph algorithms must identify when a global iteration is complete. A rank with an empty local work queue may still receive new work from another rank; local emptiness is not global termination. Termination detection therefore requires a global invariant, collective check or a carefully designed asynchronous protocol. For graph analytics, converged numerical values and absence of pending messages are separate issues. A PageRank iteration might meet a local tolerance on one rank while global residual remains large. Reduction of residual norms and compatible stopping decisions are part of correct distributed execution.

Performance measurements should include frontier distribution, visited-edge count, degree skew, communication bytes and peak memory. Graphs with equal numbers of vertices may differ in edge count by orders of magnitude. Traversing a chain has a long dependency depth even though each vertex has small degree; traversing a star may process many edges in one level but create contention around a central hub. Use several graph families, including adversarial skew, to evaluate scheduling policies. A benchmark on one synthetic random graph does not establish broad superiority.

Correctness tests use small graphs with known traversal distances, multiple paths, disconnected components, self-loops and duplicate edges. For algorithms that are meant to ignore edge multiplicity, verify the implementation actually does so. For weighted graphs, ensure the chosen algorithm supports the weight assumptions. Parallel graph processing rewards decomposition of edges and careful state ownership, but its irregularity makes hidden synchronization and load imbalance especially important.

§7.3 Machine Learning

Parallel machine learning commonly distributes data examples, model parameters or layers of computation. Data parallelism gives several workers copies of a model and different subsets of a training batch. Each worker computes gradients locally; an aggregation step combines them before the shared model is updated. Model parallelism distributes pieces of a large model across devices because a single device lacks sufficient memory or because particular operations can run concurrently. Pipeline parallelism divides model layers into stages, processing multiple microbatches in flight. These strategies may coexist, but their synchronization and statistical consequences differ.

For synchronous data-parallel gradient descent, suppose p workers hold batches of sizes n_1 through n_p. Each computes a sum of per-example gradients g_i. The correct global average gradient is the sum of all g_i divided by the total number of examples, not necessarily the unweighted average of worker-average gradients. If one worker has ten examples and another has thirty, their means must be weighted 10:30. This detail is easy to overlook when batches are uneven because the final batch of a dataset is incomplete. A numerically correct all-reduce can still produce an incorrect optimizer update if the mathematical scaling is wrong.

Communication cost can dominate when the model has many parameters but each worker processes only a small batch. If every optimizer step requires transferring a gradient vector with billions of components, network bandwidth and collective latency become central. Larger local batches reduce the frequency of gradient synchronization per processed example, but they also change optimizer behavior and potentially generalization. Gradient accumulation can simulate a larger effective batch while keeping device-memory use manageable; the learning rate and loss scaling must be consistent with the chosen aggregation. Compression of gradients can reduce bytes but introduces approximation error and extra compute, which should be evaluated on both convergence and wall-clock time to target accuracy.

Pipeline training adds scheduling constraints. If a model is divided into four stages, stage one can process microbatch two while stage two handles microbatch one. Without enough microbatches, some stages remain idle, creating a pipeline bubble. With many microbatches, utilization improves, but activation memory and scheduling complexity increase. Backpropagation dependencies complicate the schedule further because gradients flow backward through stages. Pipeline depth should therefore be chosen relative to the number of microbatches, device speed differences and allowable activation storage. A pipeline that balances nominal floating-point operations may still bottleneck on one stage's memory or communication.

Mixed precision improves throughput and reduces memory traffic on many accelerators, but model updates can become numerically unstable if small gradients underflow or large activations overflow. Loss scaling, accumulator precision and normalization schemes can mitigate risks. Parallel reductions can change floating-point order, making training runs non-bitwise reproducible even with fixed random seeds. Dropout, data shuffling and nondeterministic device kernels add further variability. Reproducibility requires recording data order, model initialization, library versions, devices and the deterministic behavior actually supported by the backend.

Asynchronous training allows workers to apply updates without waiting for every peer. This can hide stragglers but introduces stale gradients computed from older model parameters. The algorithm's convergence properties depend on delay, step size, data distribution and assumptions about the objective. It is not sufficient to claim that asynchronous updates are always faster: fewer waits can be offset by extra iterations to reach the same validation quality. Evaluate time to a defined loss or task metric, not only examples processed per second. A system can achieve extraordinary raw throughput while converging poorly.

Inference has different goals from training. Batched inference can use GPU matrix operations efficiently, but a service with strict per-request latency may not tolerate large batching queues. Dynamic batching combines requests arriving close in time, trading latency for throughput. Model replication distributes independent requests across devices, whereas tensor parallelism splits a single inference operation across devices and introduces communication. Measure tail latency, memory footprint and energy as well as average throughput. A system optimized for offline batch processing may be unsuitable for interactive use.

For verification, compare a small parallel training step against a serial reference on a fixed minibatch. Check that gradients are correctly normalized, parameter updates agree within tolerance, and no examples are dropped or counted twice. When using data loaders, ensure that distributed samplers partition examples as intended. Record validation quality after comparable compute budgets. Parallel machine learning is successful only when it accelerates progress toward the statistical task without silently changing the experimental question.

§7.4 Scientific Pipelines

Scientific computing often consists of a series of connected transformations rather than one monolithic kernel. A telescope workflow might ingest raw frames, calibrate sensor artifacts, align exposures, detect candidates, classify objects and write a catalog. Each stage has its own compute, memory, I/O and accuracy characteristics. Pipeline parallelism overlaps stages on different data batches so that, after startup, throughput can approach the capacity of the slowest stage. However, end-to-end latency for one item remains tied to the serial path through the stages unless individual stages also exploit internal parallelism.

Suppose four pipeline stages take 2, 5, 3 and 4 milliseconds per item. With one worker per stage and perfect buffering, the steady-state bottleneck is the five-millisecond stage, giving an upper throughput of about 200 items per second. The latency of one item is at least 14 milliseconds before queueing, transfers or coordination. Adding another worker to the two-millisecond stage will not improve the throughput because it is not the bottleneck. Adding a worker to the five-millisecond stage may move the bottleneck to the four-millisecond stage, assuming the slow stage can be replicated and the dataflow permits independent items. This kind of capacity reasoning is more useful than blindly parallelizing every stage.

Queueing creates another constraint. If acquisition continuously produces 250 items per second while downstream processing can sustain only 200, a queue accumulates about fifty items per second. Unless the source slows or data are discarded, memory and storage use grow without bound. Bounded queues propagate backpressure and limit resource use. In real-time observation systems, dropping low-priority data may be acceptable; in medical or scientific record processing, discarding unprocessed inputs can be unacceptable. The pipeline must state the policy explicitly. A retry policy must not accidentally process a record twice if downstream effects are not idempotent.

Data locality matters between stages. If one stage finishes a result on a GPU and the next stage also runs on the same GPU, keeping the intermediate tensor on device may avoid a costly host-memory round trip. If another stage needs CPU-only functionality, transfers are inevitable but can perhaps be overlapped with other work. Similarly, reading a large input file independently in each stage wastes I/O bandwidth. Streaming and chunked processing can reduce peak memory, but chunk boundaries must preserve algorithms that use neighboring samples or global statistics.

Workflow graphs allow tasks to execute when dependencies are satisfied. A job scheduler can track completed datasets and direct independent partitions to workers. For reproducibility, the graph should reference immutable input versions and record software, parameters and checksums. Otherwise, rerunning a failed task may read an updated input file and produce a different result while claiming to resume the same computation. Caching intermediate outputs can save expensive recomputation when only downstream parameters change, but cache keys must include all inputs that affect output. Hidden global state makes invalid caching a severe risk.

Checkpointing at stage boundaries is a tradeoff. Persisting every intermediate dataset makes recovery easy but can consume enormous storage and I/O time. Keeping only selected checkpoints reduces expense but increases recomputation after a failure. The best policy depends on stage cost, failure probability, storage price and the value of provenance. A pipeline serving regulated or high-stakes scientific conclusions may prioritize complete traceability over minimum wall-clock time. The engineering design should distinguish performance optimization from the requirements for trustworthy results.

Instrumentation should make bottlenecks visible. Measure per-stage service time, queue wait, input and output volumes, failed or retried items, worker utilization and tail latency. A high global CPU utilization is not necessarily healthy if a stage is overproducing work that piles up in a queue. Conversely, idle workers may be expected if the upstream acquisition rate is sporadic. Capacity planning follows the arrival process and required turnaround time, not merely the maximum throughput observed in a benchmark.

A scientific pipeline is successful when it produces a traceable, correct result within the user's timing and cost constraints. Parallel execution gives flexibility to overlap independent stages and distribute heavy kernels, but it also introduces queues, retries, versioning and synchronization requirements. Those concerns are part of the computational method, not separate operational afterthoughts.

THREE WORKED QUESTIONS

Step-by-Step Programming Fundamentals Solutions

Three original questions connect calculations, definitions and interpretation. Open each solution to follow the reasoning.

Solved Problem Example 7.1: Check a heat stencil

An explicit heat update uses r=0.2 and neighboring values 8 and 12 around center value 10. Compute the next value and state a stability condition for a standard one-dimensional scheme.

Form the update

u_new=10+0.2×(8−2×10+12)=10+0.2×0=10. The value remains constant in this symmetric local example.

State stability assumption

For the standard explicit one-dimensional diffusion discretization with the usual nonnegative diffusion coefficient, r≤0.5 is a familiar stability requirement.

Check parallel time levels

Every worker must read the same old time level. In-place updates may mix time levels and invalidate the numerical scheme even though the arithmetic expression looks unchanged.

Solved Problem Example 7.2: Estimate pipeline throughput

A four-stage streaming workflow needs [2,5,3,4] milliseconds per item. Compute its ideal steady-state throughput and single-item latency.

Find bottleneck service time

The slowest stage processes one item per five milliseconds, so steady-state throughput is bounded by 1/0.005=200 items per second with one worker per stage.

Find minimum item latency

One item traverses all four stages in at least 2+5+3+4=14 milliseconds without queuing, transfer or coordination.

Assess resource allocation

Adding capacity to the two-millisecond stage does not raise the 200-item/s throughput ceiling. Speeding up or replicating the five-millisecond stage can help until another stage becomes the bottleneck.

Solved Problem Example 7.3: Weight distributed gradient means

One worker has a mean gradient of 2 computed on 20 examples; another has mean gradient 5 on 30 examples. Compute the global average.

Recover local sums

The workers contribute sums 20×2=40 and 30×5=150. Total gradient sum is 190 over 50 examples.

Compute result

The correctly weighted global mean is 190/50=3.8; the unweighted mean 3.5 is inaccurate because the shard sizes differ.

Explain application

All-reduce sums and counts, then divide. Scaling mistakes change optimizer steps and can make comparisons across worker counts scientifically invalid.