gpumetropolis 0.7: from a fast sampler to a full joint model
The 0.4 guided tour left gpumetropolis as a fast random-walk sampler with a decision layer and an honest benchmark. The releases since then added two things the tour did not have: an engine and a destination. The engine is the gradient, exposed by differentiating the compiled log-density, which turns on a Langevin sampler and a conjugate exact path. The destination is the full joint distribution, reached in 0.7 by pairing automatic marginal selection with the copula layer of 0.6. This post is the map of that arc, release by release, ending where 0.7 lands: you hand it a data frame and it models the joint, marginals and dependence, each chosen by the data.
The engine: a gradient from the bytecode
The package compiles each model’s log-density to its own bytecode. Version 0.5.2 differentiates that bytecode in reverse mode, so the gradient of the compiled density is exact, validated against finite differences to 1e-10 relative error, and it is JIT-compiled to straight-line native code next to the density itself. A gradient step therefore costs a small multiple of a plain evaluation rather than a separate numerical pass.
With an exact gradient in hand, method = "mala" becomes available: the Metropolis-adjusted Langevin algorithm, whose proposal drifts along the gradient of the log-posterior and whose acceptance carries the asymmetric-proposal correction. The drift is truncated at a whitened norm of 2 sqrt(d) (ROBERTS; TWEEDIE, 1996), so a chain started far from the mode walks in instead of freezing, and because the correction uses the truncated mean the invariant distribution stays exact. The warmup targets the Langevin acceptance optimum of 0.574 (ROBERTS; ROSENTHAL, 1998). On the high-dimensional case the earlier random walk had conceded, a logistic regression in twenty-one parameters, the Langevin path reaches about 30,000 effective draws per second against roughly 10,000 for a gradient sampler with a compile step and about 3,000 for the plain random walk; on the applied cases it multiplies the earlier wins. The same 0.5 line also learned a full-covariance adaptive random walk, which pools the covariance estimate across chains rather than trusting many noisy per-chain ones, and a conjugate exact fast path (gpum_lm) that recognises a Gaussian linear model under a Normal-inverse-Gamma prior and draws the closed-form joint posterior directly, no chain at all, at millions of effective draws per second where a Gibbs specialist reaches a fraction of that.
Beyond independent rows
The 0.5.3 release states the principle that lets the same per-row engine handle dependent data: the engine’s per-row sum is the exact log-likelihood whenever each row is the conditional density of its observation given the ones before it. That covers independent non-identically distributed rows, which every regression already is, and stationary Markov dependence. gpum_ts_model() assembles the lagged conditional design of a Markov model, validated against arima() on simulated AR(1) and on the Nile series, and gpum_lfo() does exact leave-future-out cross-validation with an expanding window, the comparison tool that respects temporal order where the pointwise exchangeability behind WAIC and PSIS-LOO fails (BURKNER; GABRY; VEHTARI, 2020). The boundary is stated in the same place: latent recursions such as GARCH and state-space models, and matrix-coupled likelihoods such as Gaussian processes, are outside the per-row language, and the documentation names the right tools for them rather than leaving the limit quiet.
The destination: modelling the joint distribution
Version 0.6 added the copula. gpum_copula() fits a copula to the pseudo-observations of two columns, average ranks rescaled to the open unit square, so the fit is invariant to the marginal shapes and isolates the dependence structure. Four families cover the qualitative range: Gumbel for upper-tail dependence, Clayton for lower-tail, Frank for symmetric tail-independence, and Gaussian for the elliptical case that also carries the sign of the dependence. The family is a hypothesis about the tails, so family = "auto" fits every candidate and returns the one preferred by predictive comparison, withholding any whose chains do not converge.
Version 0.7 supplies the other half. gpum_fit_catalog() takes a numeric column, detects its support and modality, fits every eligible parametric family with the sampler, and ranks them by predictive comparison. Ten families ship: normal, Student-t, logistic and Laplace on the real line; gamma, lognormal, Weibull and exponential on the positive half-line; beta on the unit interval; and a two-component Gaussian mixture for a bimodal column. Every family is written in an unconstrained parametrisation, the log of each positive parameter, so the sampler never proposes an invalid scale or shape, and every log-density carries its normalising constant so the comparison is on a common scale. The winner’s fitted cumulative distribution is exactly the probability-integral transform the copula needs, returned by marginal_cdf().
Put the two together and you have both halves of Sklar’s theorem (SKLAR, 1959): the explicit marginals from the catalogue and the copula that binds them model the full joint distribution. The 0.7 case-study vignette walks it end to end on the airquality pair, each column given its own selected marginal before the copula couples them, and on the Old Faithful waiting times, where the catalogue selects the mixture and recovers the two modes. The bimodal fit is the delicate one: the two components start at the data quartiles with their means held by a data-anchored prior, which removes the empty-component degeneracy that otherwise sends a mean to infinity, and the draws are relabeled so the lower-location component is always reported first.
One log-gamma, three backends
Five of the ten families need the log-gamma of a shape parameter in their density, so 0.7 adds lgamma as a new operation across the whole stack: the interpreter, the cranelift JIT, the reverse-mode automatic differentiation with its derivative the digamma, and the GPU kernel. The point of the package has always been that a single portable formula runs on CPU, CUDA and Vulkan, so a new operation is only real when it produces the same inference on all three. On the CPU path the value matches R to 1e-15 and the gradient to 1e-9. On the GPU the kernel is a Stirling approximation in single precision, runtime-verified on an RTX 5090 on both engines, the CUDA backend and the Vulkan backend through wgpu, each matching the CPU backend to Monte-Carlo noise on the normal, gamma and beta models, with a mean absolute difference below 1.6e-4. Where the density is expensive and the chains are many, the accelerator earns its place: CUDA delivers an eight-fold speedup at 2,048 chains on a gamma target with eight thousand observations.
That same hardware exposed a failure worth catching well. A CUDA kernel cannot compile for a GPU newer than its toolkit, and the raw runtime error is an opaque nvrtc: invalid value for --gpu-architecture string. Version 0.7 catches the launch failure and follows it with the probable cause and the fix: the toolkit is older than the GPU requires, so rebuild against a matching toolkit, or run with backend = "vulkan" or backend = "cpu". An error message that names the cause is part of the portability claim, not a nicety beside it.
Where it sits
The arc from 0.4 to 0.7 is a single move: give the fast sampler a gradient, then point it at the joint distribution instead of a single parameter vector. You can now hand gpumetropolis a data frame and get back a model of the joint, each marginal chosen from ten families by predictive comparison and the dependence chosen from four copulas the same way, computed by one portable formula that runs identically on a laptop CPU and a Blackwell GPU. The boundaries remain named rather than hidden: latent-recursion and matrix-coupled likelihoods are still outside the per-row language, and the roadmap, not this post, is where they live.
The package is on R-universe with the vignettes that work each of these cases in full, and the source is on GitHub: github.com/pcbrom/gpumetropolis. I write about this kind of work, scientific method, statistics, and AI applied with rigor, on LinkedIn: linkedin.com/in/pcbrom.
References
BURKNER, P.-C.; GABRY, J.; VEHTARI, A. Approximate leave-future-out cross-validation for Bayesian time series models. Journal of Statistical Computation and Simulation, v. 90, n. 14, p. 2499-2523, 2020. DOI: 10.1080/00949655.2020.1783262.
ROBERTS, G. O.; ROSENTHAL, J. S. Optimal scaling of discrete approximations to Langevin diffusions. Journal of the Royal Statistical Society: Series B, v. 60, n. 1, p. 255-268, 1998. DOI: 10.1111/1467-9868.00123.
ROBERTS, G. O.; TWEEDIE, R. L. Exponential convergence of Langevin distributions and their discrete approximations. Bernoulli, v. 2, n. 4, p. 341-363, 1996. DOI: 10.2307/3318418.
SKLAR, A. Fonctions de repartition a n dimensions et leurs marges. Publications de l’Institut de Statistique de l’Universite de Paris, v. 8, p. 229-231, 1959.