GPU support

Every process carries its sampling backend as a type parameter and is re-typed with .on::<B>() (see backends). This page is the inventory: which processes have a device kernel today, on which backends, and what the rest do.

The support matrix

Every process reaches a device: Fgn and Fbs through pipelines of their own, the other 130 through the engine. The 130 processes on the engine cover the diffusions, the stochastic-volatility models under both their Euler and quadratic-exponential schemes (Heston, HestonLog, DoubleHeston, Heston2D, Sabr, Bergomi, HestonStochCorr, BatesSvj, Hkde), the short-rate models including the curve-driven ones (HullWhite, HullWhite2F, HoLee, BlackKarasinski, CirPlusPlus, Cir2F, DuffieKan, DuffieKanJumpExp, Vasicek, FVasicek), the fractional processes, the jump models (MjdLog, BatesSvj, Hkde, HawkesJD, DuffieKanJumpExp, and the ones generic over a user-supplied size law — Merton, Kou, LevyDiffusion, Bates1996, JumpFou, JumpFOUCustom — whose ScalarNormal or ScalarExp law the engine recognises at runtime and draws in the kernel, Bates's as the product ∏(1 + y) − 1 its multiplicative jump needs; any other law keeps the process on the host through the same calls), the Lévy processes whose increment is one draw with no rejection (Ig, Nig, IGSubordinator, PoissonSubordinator, AlphaStableSubordinator, GammaSubordinator, TemperedStableSubordinator, Vg, BilateralGamma, BilateralGammaMotion, Bns), the discrete-time conditional-variance models at order (1,1) (Arch, Garch, GjrGarch, Agarch, Egarch), the first-order lag models (ARp, MAq), the noise generators (Wn, Gn, Cgns, Cfgns), the correlated fractional pair (Cfbms, Cfou), the bounded correlation processes (TengSCP, VanEmmerich, TransformedOU), the two models whose coefficients travel as curves rather than parameters (MultifactorSabr's term structure, Hjm's six time-dependent coefficients), the matrix-valued models whose rows are independent copies of one small family (Bgm's forward rates, Adg's affine factors, WuZhangD's rate–variance pairs — one launch per row), the Cholesky-correlated baskets of up to four streams (MultiGbm, Mcgns — one launch, the factor folded into the family's ten lower-triangle slots; a wider basket stays on the host whatever the backend), the one- and two-factor MultifactorHeston (which rides the double-Heston family with a zeroed second factor when it has one; three or more stay on the host), the regime-switching diffusion with up to four regimes (the one-step transition matrix is exponentiated once on the host and travels as cumulative thresholds; five or more regimes stay on the host), the Riemann-Liouville family under the Markov lift of its rough kernel (RlFBm, RlFOU, RlBlackScholes, RlHeston, and Volterra's lift branch — the lift's node states live in the kernel's own per-path arrays, see below; Volterra's reference-convolution branch rides the history block instead, the kernel tabulated on the grid as its weights), the two lift-stepped Volterra models (VolterraSquareRoot, and GaussianPolynomialVolatility up to eight coefficients) and the event-indexed processes in their count mode (Poisson, CompoundPoisson, and CustomJt / CompoundCustom when their inter-arrival law is the scalar exponential), one arrival per step, and the models whose step reads their whole past through the history block — the rough-volatility models with a memory integral (RoughHeston, FBatesSvj, RoughBergomi), the linear fractional stable motion Lfsm, and the ARIMA and SARIMA recursions at any order, run as the moving-average filter of their innovations — each on the host's own scheme, for grids up to 512 points; and the tempered-stable processes built from their shot-noise series (Cgmy, Cts, KoBoL, Rdts), whose terms the series block draws per path and sums into the grid cells they fall in, the host's sort replaced by addressing, again up to 512 points; and Hawkes in its count mode, one event per step by the exact two-uniform recursion of Dassios and Zhao, the host's thinning loop replaced by a closed form; Ctrw under exponential waits, whose arrivals per cell are the Poisson count the kernels draw, with normal, Rademacher or symmetric-stable jumps; and Lmm with up to four forwards in one launch, the spot measure's drift coupling and the reset freezes written out over the four slots; and InverseAlphaStableSubordinator, whose table block builds the direct stable subordinator per path on its own grid — doubling the extent until it reaches the horizon, as the host does — and reads the inverse clock off it by interpolation at every step; and Svcgmy, CGMY under a CIR variance, whose variance takes the exact square-root step — a non-central χ² as the square of a shifted normal plus the frame's gamma draw — and whose series terms are sized in the step of their own cell against the variance there, the one series clause written live; and the two-dimensional Wishart process by its exact Ahdida–Alfonsi step — the three maps of the law identity folded on the host, two coordinate steps each a shifted normal's square plus the frame's central χ² draw and a Gaussian move of the off-diagonal entry — at a degree of three or more, where the path stays inside the cone and the host's rank-adaptive branch is never taken; and MultivariateHawkes in count mode with two components and one decay per target, as the superposition of each target's exact excess clock (Dassios–Zhao) with the baselines' joint Poisson clock, the event being the earliest of the three — the two extra uniforms come from the noise components, half the sum of two squared normals being a unit exponential and one minus its negative exponential a uniform. The two-dimensional field Fbs is the one other device process off the engine: its sample is a transform, not a recursion, and it runs through a pipeline of its own behind SheetBackend<T>, described below.

What stays on the host, and why. The engine steps a bounded forward recursion over four scalar state slots, so what rules a process out is unbounded work per step, unbounded state, or an output that is not a fixed grid:

• A history integral over every earlier step, on a grid longer than 512 points (Volterra's reference branch and any history family above that length): the step sums (t_i - t_j)^(H - 1/2) z_j over all j < i, which is work that grows with the step index. The Markov lift is the answer where it applies — the whole Rl* family and Volterra's rough kernels ride it — and the history block runs the direct convolution itself, exactly as the host does, for grids up to 512 points; that is how RoughHeston, FBatesSvj, RoughBergomi and Lfsm reach a device. Longer grids stay on the host, which the process's device_ready() decides.
• A variance with less than one degree of freedom (Svcgmy at 4κη/ζ² < 1): the exact square-root step then draws a Poisson mixture the kernels do not carry, and the process stays on the host; at one degree or more the step is the square of a shifted normal plus a central χ², which the frame's gamma draw supplies, and the series block sizes the terms live.
• A runtime dimension that couples its rows, above the slots (Wishart above d = 2, Lmm / MultiGbm / Mcgns above four): the rows share a drift or a correlated noise, so they cannot be stepped as independent launches the way Bgm, Adg and WuZhangD are, and a launch carries four slots and four shocks; within them they ride one launch, above they stay on the host whatever the backend. Wishart also stays below a degree of three, where its exact step may take the rank-adaptive branch — or, below two, the Poisson mixture — that a family does not express.
• More time-varying coefficients than the eight the kernel binds: none today. MultifactorSabr needed three, Adg five and Hjm six, and all three now declare them — the first turned a per-step bucket search into a buffer read, the last seven time-only callables into six tabulated curves.
• A coefficient written as a Rust closure or a Python callable (Cheyette's sigma(t, x), VolterraSde's drift and diffusion, HestonSlv's leverage): a closure of time alone is tabulated as a curve the way Cir2F's phi(t) is, but one that reads the state has no representation a kernel can evaluate, and neither has a coefficient tabulated as a Grid2D — a calibrated leverage surface, say. Written as an Expr instead — Expr::lit(0.01) + Expr::x().abs() * 0.4 — the same coefficient compiles to a short program the kernel interprets, and the process runs on the device; see "Coefficients as expressions" below.
• Point processes on their own clock (Hawkes and MultivariateHawkes in their horizon modes, and MultivariateHawkes above two components or with per-pair decays): the output length is the number of events, not the grid. Hawkes in count mode is one event per step and reaches the engine through the exact exponential-excitation recursion, and MultivariateHawkes in count mode with at most two components and one decay per target by superposing each target's exact excess clock with the baselines' joint Poisson clock — per-pair decays make a sum of exponentials at different rates whose first arrival has no closed form, and a third component would need a fifth state slot, so both stay on the host; Ctrw is a grid process and reaches it under exponential waits, while its gamma, inverse-Gaussian and stable waits have no per-cell Poisson count and stay on the host. Poisson, CompoundPoisson, CustomJt and CompoundCustom split here — their count mode is a fixed-length running sum and reaches the engine, their horizon mode does not.
• A per-step search over a table finer than 512 points (InverseAlphaStableSubordinator with more u_steps than that): the table block keeps the direct subordinator in a per-path array of 512 values and scans it; a finer table stays on the host.
• More regimes than the four the chain family carries (RegimeSwitchingDiffusion at five or more): the family holds four volatilities and four rows of thresholds; larger chains stay on the host whatever the backend.
• A lag order above one (ARp, MAq or the GARCH family built above order one): the recursion reads back as many past values as the order and a step holds four slots. At order one they reach a device; above it the process stays on the host, as its device_ready() reports. Arima and Sarima are the exception: a linear recursion of white noise is a moving-average filter, so the history block runs them at any order as the convolution of the innovations with the recursion's own impulse response.
• A sheet whose embedding sides are not powers of two (Fbs with m − 1 or n − 1 not a power of two): the sheet pipeline's radix-2 transforms take embedding sides 2(m − 1) and 2(n − 1) that are powers of two; any other grid stays on the host whatever the backend.
• More variance factors than a launch's four shocks carry (MultifactorHeston at K >= 3): each factor takes two shocks, so one and two factors ride the double-Heston family and three or more stay on the host whatever the backend.

Coefficients as expressions

Cheyette, VolterraSde and HestonSlv take their coefficients as [Fn2D] callables of time and state. A Rust closure or a Python callable runs on the host only — nothing can evaluate it in a kernel — and the device build keeps such a process on the host, bit-identically to the Cpu build. Written as an [Expr] the same coefficient reaches the device:

use stochastic_rs::traits::Expr;

// σ(t, x) = 0.01 + 0.4 |x|, a displaced local volatility
let sigma = Expr::lit(0.01) + Expr::x().abs() * 0.4;
let cheyette = Cheyette::<f32, _>::new(f0, 0.5, sigma, 252, Some(1.0), seed).on::<Cuda>();

Expr::t() and Expr::x() are the time and the state, Expr::lit a constant; the arithmetic operators, sqrt, exp, ln, abs, tanh, powf, max and min build the tree — the vocabulary the engine's own function table has, so every node means one thing on the host and in every kernel. From<Expr> for Fn2D compiles the tree to a postfix [Program] of at most 62 operations on a stack of at most 8, checked at compile time on the host. The host evaluates that program wherever it would have called the closure; a launch uploads it — two programs at most, [len₁, len₂, (opcode, constant)…] in one small buffer — and before every step the kernel runs each program at the step's time and first state slot, handing the values to the family as pv and pv2. The CheyetteLocalVol family reads its volatility as pv, the HestonSlv family its leverage; the VolterraProgram family, which steps under the kernel's Markov lift, reads its drift and diffusion as pv and pv2. Each kernel — the C body CUDA and Metal render — carries the same sixteen-opcode interpreter, and the parity suite runs a probe program with every opcode in use through all of them against the host's Program::eval.

What a family is

A process reaches a device by declaring which of the engine's 121 families it steps, not by writing a kernel. A family is one entry in stochastic-rs-stochastic/src/euler/families.rs — its parameters, its state components, its noise components, one step expression per component and one report expression — and the same tokens generate the host step, the C body the CUDA and Metal kernels render. Adding a process is a declaration plus an EulerCoefficients (one component) or EulerSystem (more than one) impl; no kernel is edited.

The engine's kernels carry four state slots, four independent noise components, twenty scalar parameters, two uniforms, a Poisson jump count, the sum of that step's jump sizes — under a normal, double-exponential, tempered-stable, Rademacher or symmetric-stable law, the first two also compounded as ∏(1 + y) − 1 for a multiplicative jump or taken one per step for an event-indexed path — and up to two Gamma draws. A family with a lift clause runs under the Markov lift of a rough Volterra kernel: the frame keeps up to 176 nodes' two state vectors per path, evaluates the family's drift, diffusion and shock before the step and hands it the lifted value as lv. A family with a history clause keeps, per path, the values it pushes each step — up to 512 grid points — and reads their convolution with a weight curve the launch tabulates, indexed by lag, as cv: the exact O(n²) scheme of a moving average over the whole past, which is a Volterra integral on the grid, an ARMA filter or a fractional stable motion depending on the weights. A family with a series clause has the frame draw, per path and before the steps, as many shot-noise terms as the launch names — a unit-rate arrival, an exponential, two uniforms and an arrival time each — size them with the family's own expression and sum them into the grid cells their times fall in, up to 512; the step reads its cell as sj. A family with a table clause has the frame build, per path, a monotone table of up to 512 values from its increment expression, doubling the extent until the table reaches the horizon, and hands each step the interpolated inverse at its time as iv. A launch may carry up to two programs — coefficients written as Expr — which the frame evaluates before every step at the step's ct and first state slot and hands to the family as pv and pv2. A launch binds up to eight time-varying coefficients, which a family names ct and ct1 through ct7; it pays one buffer read per declared curve per step, so a model with one still pays for one. The first two noise components may come from a fractional pipeline instead of the kernel's own hash stream: one embedding produces both when the rows share a Hurst exponent, which is what a correlated fractional pair needs. A process whose first grid point is itself a draw — a conditional-variance model's series, a noise generator — asks the launch to step before writing it. A model that wants correlated shocks draws the components independently and correlates them in its own step, which is what the host samplers do.

Feature flags: metal (Metal), cuda (Cuda), accelerate (Accelerate). A marker only exists when its feature is compiled, so selecting a backend a process does not support — or one the build does not carry — is a compile error, never a silent fallback. The same holds for precision: the capability traits take the scalar (FgnBackend<T>, SheetBackend<T>, EulerBackend<T>) and a device implements them only for the precision its kernels compute in, so Gbm<f64> or Fgn<f64> on Metal does not compile. Nothing is computed in f32 behind an f64 type.

What the device path does

• fGN noise. The circulant-embedding FFT that produces fractional Gaussian noise runs on the device, one launch per sample_par batch. Fgn returns it directly.
• Programs. A coefficient written as an Expr travels to the kernel as postfix code and is interpreted before every step — Cheyette's local volatility, VolterraSde's drift and diffusion, HestonSlv's leverage — so the device runs the very same coefficient the host does, not a tabulation of it.
• Sheets. Fbs runs its two-dimensional circulant embedding as a pipeline beside the fGN one. The embedding's eigenvalue roots are computed once on the host, when the process is built; per batch the device draws the complex Gaussian noise, scales it, transforms the rows, transposes, transforms the columns and reads out the leading m × n block less its corner plus Stein's linear correction √(2c₂)(t₁Z₁ + t₂Z₂). The radix-2 transforms take embedding sides that are powers of two, so a grid with m − 1 and n − 1 powers of two runs on the device and any other stays on the host, bit-identically to the Cpu build. A batch produced in chunks equals one launch sheet for sheet, as with the fGN pipeline.
• Fractional diffusions. Fbm, Fou, Fcir, Fgbm, FJacobi, FVasicek, Cfgns, Cfbms and Cfou hand those increments to the Euler engine, which runs the whole recursion in the kernel. They declare no kernel of their own: Fou is the same Ornstein–Uhlenbeck family a Gaussian Ou uses, Fgbm the same geometric one, and a step that multiplies by an increment does not care whether the increment was hashed or produced by an FFT. The two kernels meet in device memory: the pipeline hands its own buffer to the Euler launch instead of reading the increments back, so a batch crosses the bus once, for the paths. Measured on an M-series Metal device against the host: fOU with 2000 paths of 512 steps went from 101 ms to 5.8 ms, fBM with 4000 paths of 1024 steps from 252 ms to 19.2 ms. Keeping the increments on the device is most of that: through host memory the same runs took 25 ms and 93 ms.
• Euler engine. The whole Euler–Maruyama recursion of a path runs in one kernel thread; sample_par(m) is one launch for all m paths, sample() launches a single path. A system writes one plane per component, so a Heston draw comes back as the two arrays it always did. CIR uses full truncation (Lord, Koekkoek & van Dijk 2010) in the kernel. The CUDA and Metal kernels are one C body rendered per language, so Metal and Cuda agree with each other seed for seed up to f32 libm rounding; euler::family_parity launches every declared family on the device present and holds it to the host step.
• Host backends. Cpu is the process's own sampler; Accelerate is a CPU device that only differs for the fGN family (vDSP FFT). For every other process both are the same code path.

Large batches

A device handle holds at most batch_budget bytes of path data per launch (1 GiB by default, STOCHASTIC_RS_DEVICE_BATCH_BYTES, or Cuda::default().with_batch_budget(n)). A larger sample_par(m) runs as several launches whose union is bit-identical to one launch, because the kernels hash the global path index; sample_map(m, f) maps each chunk in parallel on the host before the next chunk is launched, so the batch never has to fit in memory at once. On native CUDA the chunks alternate between two streams: one copies back through pinned memory while the next computes. The fGN devices chunk the same way, one seed per batch and an element offset per chunk, and keep the FFT plans and buffers of the last four (n, m) shapes, so alternating two sizes does not re-plan on every call.

More than one GPU

The handle names the device: Cuda::new(1) (CUDA ordinal), Metal::new(1) (index into Device::all()); process.with_backend(handle) puts a process on that handle, .on::<B>() on the backend's default handle, and handle.probe() reports the device opened. The defaults read STOCHASTIC_RS_DEVICE and STOCHASTIC_RS_DEVICE_BATCH_BYTES. From Python the name carries the ordinal: device="cuda:1", probe_device("metal:0").

When the device is missing

handle.probe() on any backend handle returns Ok(DeviceInfo) with the device's name and precisions, or Err(DeviceError) saying why it cannot be used (Unavailable, Compile, Launch). Every process has try_sample() and try_sample_par(m) on ProcessExt, returning the same error instead of panicking (always Ok on the host devices, bit-identical to sample*); the plain sample* calls panic with the error's message and never fall back to the host. From Python a device failure is a RuntimeError.

A configuration the kernels cannot carry is a different matter and follows one rule: the process samples on the host, bit-identically to the Cpu build, whatever backend it is on, and device_ready() — a ProcessExt method, true by default — says so beforehand. Nothing panics and nothing changes law; the "What stays on the host" list above is exactly the set of configurations that answer false.

Precision and reproducibility

The device seed is drawn from the process's own seed source, so the device follows the host's rule exactly: two processes built from the same Deterministic seed value give the same paths, consecutive calls on one process advance the stream into independent paths, and an Unseeded process draws fresh entropy every call. Host and device paths agree in distribution, not bit for bit, because the device draws its own stream.

From Python

Every class whose process reaches a device takes device= (PyGbm(..., dtype="f32", device="metal")), checked and probed at construction. A class also answers device_ready() and device_fallback() for the configuration it holds, so a device= that the kernels cannot honour is visible rather than merely slow. PyGbm(0.05, 0.2, 253, x0=100.0, t=1.0, seed=7, device="metal", dtype="f32").sample_par(m) runs the Euler engine with device one of "cpu", "accelerate", "cuda", "metal" — every name is one concrete backend, there is no "pick whatever is compiled" name; one the build does not carry raises ValueError with a rebuild hint. The array is float64 from the CPU and native CUDA paths and float32 from Metal, whose kernels compute in single precision. The published wheels are CPU-only; the device backends come with a source build — maturin develop --features metal on a Mac, --features cuda on an NVIDIA machine.

Trying the CUDA back-end without a CUDA machine

notebooks/colab_cuda_check.ipynb opens in Google Colab: on a T4 runtime it installs Rust, clones the repository, runs the native CUDA tests (the fGN pipeline and the Euler engine), builds the Python wheel with device="cuda" and plots what the GPU produced. The Rust tests build in 10–20 minutes on Colab's two cores, the wheel in another ten; the tests and the plots take seconds.

A run on a T4

notebooks/runs/colab_cuda_check_T4_2026-09-03.ipynb is the executed copy of a run on 2026-09-03 (Tesla T4, driver 580.82 / CUDA 13.0, toolkit 12.8; it predates the rename of the backend, so its cells, its output and its handle names still read cuda-native / CudaNative; GitHub's own viewer may not render Colab outputs, nbviewer does). All 12 native-CUDA tests passed, including the two that assert a chunked batch equals one launch path for path, in f32 and f64; the wheel built for CPython 3.13 and its device tests passed; probe_device reported {'backend': 'CudaNative', 'name': 'Tesla T4', 'precisions': ['f32', 'f64']} (the handle is Cuda and reports "Cuda" since the rename).

The variance of fractional Brownian motion sampled on the GPU sits on the t^{2H} lines for H = 0.3, 0.5, 0.8 exactly where the CPU samples sit, and the Euler engine's GBM terminal mean over 100k paths was 105.15 against the 105.13 forward. Wall time per call, warm, against Colab's two CPU cores:

The Euler case is slower on the T4 in f64 because that card executes double precision at 1/32 of its single-precision rate, and the 200 MB result is copied back and re-laid as rows on the host; dtype="f32" is the setting for the Euler engine on inference-class GPUs. The fGN pipeline is FFT-bound and wins already at these sizes; the gap widens with m.

Adding a kernel

A process joins the Euler engine inside this crate: one family entry in euler/families.rs — parameters, state, noise, step and report expressions, plus a lift, history, series or table clause where the recursion needs one — generates the host step, the C body the CUDA and Metal kernels render and the process declares that family through the engine's launch hooks and switches its backend_switch! line to via euler. The hooks, the family specification and the slot caps are the engine's own machinery: public so the crate's tests and probes reach them, hidden from the rendered documentation and outside the stability promise. What is public and stable is the capability a bound names, EulerBackend<T>, together with .on::<B>(), with_backend(handle), device_ready(), probe() and the try_* calls. Widening a process's bound from HostBackend to EulerBackend is not a breaking change for its callers.

A coefficient a process takes as a closure has no representation a kernel can evaluate; written as an Expr it compiles to a program the kernel interprets, which is how Cheyette, VolterraSde and HestonSlv run on a device. Sde, which drives its own loop over caller-supplied closures, is the one process whose .on::<B>() re-types without changing where it samples.