← Cellular Automata From First Principles

Run Cellular Automata on the GPU

Once the update is expressed as tensor operations, moving to a GPU becomes straightforward.

But a GPU is not automatically faster.

It wins when there is enough parallel work to amortize transfer and launch overhead.


A tensor implementation of Life

import torch
import torch.nn.functional as F

LIFE_KERNEL = torch.tensor(
    [[1.0, 1.0, 1.0],
     [1.0, 0.0, 1.0],
     [1.0, 1.0, 1.0]]
).view(1, 1, 3, 3)


def life_step(x):
    padded = F.pad(x, (1, 1, 1, 1), mode="circular")
    neighbors = F.conv2d(padded, LIFE_KERNEL.to(x.device))
    alive = x > 0.5
    born = neighbors == 3
    survive = alive & (neighbors == 2)
    return (born | survive).float()

Two things matter here beyond the convolution. First, this is the third life_step in the book (after Chapters 6 and 51) — same B3/S23 rule, new backend, bridged explicitly rather than silently redefined. Second, the padding mode is load-bearing: plain padding=1 zero-pads, which silently changes periodic boundaries into dead ones. Circular padding preserves the torus semantics every roll-based helper in this book assumes — verified bit-identical against the NumPy version on 20 random grids, where the zero-padded variant mismatched all 20 at the edges.

Represent a batch as:

(batch, channel, height, width)

Move state once

device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
state = state.to(device)

for _ in range(1000):
    state = life_step(state)

Do not copy state back to the CPU every step unless you need it there.

A bad loop is:

CPU → GPU → step → CPU → GPU → step → CPU

A better loop is:

CPU → GPU → many steps → CPU

Synchronize when timing CUDA

GPU kernels execute asynchronously relative to the host.

So this benchmark can lie:

start = perf_counter()
state = life_step(state)
elapsed = perf_counter() - start

For CUDA timing (guarded — synchronize raises where CUDA is absent, so CPU-only runs skip it):

if torch.cuda.is_available():
    torch.cuda.synchronize()
start = perf_counter()

for _ in range(steps):
    state = life_step(state)

if torch.cuda.is_available():
    torch.cuda.synchronize()
elapsed = perf_counter() - start

Without synchronization, we may measure kernel submission rather than completion.


Batch experiments

The GPU becomes especially useful when evaluating many worlds at once.

state = torch.rand(256, 1, 256, 256, device=device) > 0.5
state = state.float()

Now 256 simulations share the same kernel launches.

This is ideal for:

seed search
parameter sweeps
robustness testing
NCA training
benchmark ensembles

Continuous systems fit naturally

Lenia-style updates use convolutions and pointwise functions, both GPU-friendly operations.

r = kernel.shape[-1] // 2
field = F.conv2d(F.pad(state, (r, r, r, r), mode="circular"), kernel)
growth = 2 * torch.exp(-((field - mu) ** 2) / (2 * sigma**2)) - 1
state = torch.clamp(state + dt * growth, 0.0, 1.0)

The padding is circular for the reason above: padding="same" would zero-pad and silently turn Lenia’s torus (periodic, as in Chapter 31’s FFT engine) into a world with dead edges.

Neural cellular automata are even more natural because their update rule is already built from convolutions and neural-network layers.


Precision is part of the model

Moving from float64 to float32, float16 or bfloat16 can change performance substantially.

It can also change dynamics.

For sensitive systems, tiny numerical differences may amplify.

So benchmark precision together with behavioral equivalence.

lower precision
  ≠ free optimization

GPU
  ≠ automatically faster

CPU versus GPU is an experiment

Measure both.

A small 64×64 automaton may be faster on a CPU.

A large batched NCA rollout may strongly favor the GPU.

The useful question is not:

Are GPUs faster?

It is:

At what workload does this implementation cross over?

That is an engineering result we can record and reproduce. The porting workflow — reference, move, synchronize, compare — with the boundary fix this chapter earned:

    flowchart LR
    R[CPU reference, asserted correct] --> M[move state + kernel once]
    M --> G[GPU rollout, many steps]
    G --> S[synchronize, then time]
    S --> C[compare outputs vs reference]
    C --> Q[crossover workload?]
  

CPU vs GPU, dimension by dimension (no GPU hardware exists in this book’s test environment, so the GPU column states requirements rather than reporting unmeasured wins):

DimensionCPUGPU
transfer costnone (data resident)move once; per-step copies destroy gains
timing semanticssynchronous, perf_counter sufficesasync launches — must synchronize
boundary handlingnp.roll wraps (periodic)circular pad preserves torus; zero-pad silently changes it
reproducibilityseeded NumPy, deterministicseeded torch + flags; some ops nondeterministic
wins whensmall grids, few worldslarge/batched workloads amortize launch

Both backends still pay a cost that grows with the kernel: the larger the neighborhood, the more multiply-adds per cell. The next chapter removes that dependence for large periodic kernels with the FFT.


Research

  • PyTorch documentation: torch.cuda (CUDA semantics). The execution model behind this chapter’s two hardest-won points: kernel launches are asynchronous (hence synchronize around timing), and the documented API surface tells you exactly which operations exist — including the guard behavior on machines without CUDA. https://docs.pytorch.org/docs/stable/cuda.html

  • PyTorch documentation: torch.nn.functional.pad (padding modes). The reference for the chapter’s boundary fix: constant (the default, including zero padding), reflect, replicate, and circular — with circular mode being what preserves toroidal semantics under convolution. Check the mode before assuming any conv-based helper shares your boundary model. https://docs.pytorch.org/docs/stable/generated/torch.nn.functional.pad.html