Same Workload, Opposite Throughput: A Minecraft-Shaped Speedup for Scientific Interpolation

A 57-hour problem

The cosmic density reconstruction study needed an adversarial validation campaign: many full reconstructions of a 128³ voxel grid, across many random seeds, to make sure the headline results weren't a fluke of one draw. The bottleneck was the radial basis function (RBF) interpolation step that reconstructs the grid from sparse observations.

The naive pipeline, Python with scipy's RBFInterpolator, was correct but slow. Each of roughly 2 million unobserved voxels got its own treatment: find the k=500 nearest observed neighbors with a KD-tree lookup, build a local RBF system, solve it. That adds up to two million independent tree queries and two million independent linear-algebra solves. On a 32-core workstation, one full reconstruction took about 4.3 hours, and the planned validation campaign penciled out to an estimated 57 hours of wall clock. The eventual four-experiment adversarial campaign, 380 reconstructions in total, would have needed on the order of 1,600 machine-hours at that rate. That is not a realistic wait on consumer hardware.

The recognition

The restructuring idea didn't come from the numerical-methods literature. It came from Minecraft.

Voxel game engines never process blocks individually. The world is divided into spatial chunks (16³ blocks in Minecraft's case), and each chunk shares its spatial data structures, such as lighting maps, adjacency tables, and collision meshes, across every block inside it. The engine processes chunks and reuses shared context across spatially adjacent elements, because per-block processing would be hopeless.

The interpolation workload has the same shape. Voxels that are spatially close share most of their k-nearest neighbors: the lookup for voxel (x, y, z) returns nearly the same neighbor set as (x+1, y, z). Chunk the domain and do one KD-tree lookup per chunk instead of one per voxel, and the number of independent lookups drops from about 2 million to about 4,096. Once every voxel in a chunk shares the same neighbor set, all of the chunk's local RBF systems can be batched into a single GPU linear-algebra call: one solve with hundreds of right-hand sides instead of hundreds of separate solves.

This trick has a name

This is the part that turned a research paper into a blog post.

When I ran this write-up through adversarial prior-art review, the verdict was blunt: spatially binning scattered data into a uniform grid to expose tile-local GPU locality is thoroughly established practice. NVIDIA's canonical particle-simulation sample documented the pattern in 2008 (Green 2008). It was applied to molecular modeling the same year (Rodrigues et al. 2008; Anderson et al. 2008), refined for electrostatics shortly after (Hardy et al. 2009), and codified as the named "Input Binning" pattern in Kirk and Hwu's standard GPU textbook (2011), which notes that "small bins also serve as tiles for locality." It has since become standard in GPU molecular dynamics, SPH fluids, MRI gridding, and nonuniform FFTs.

To be clear about what this post is and isn't: the technique is not new, and I'm not claiming it is. What I have is an application story. Scattered-data interpolation for cosmological reconstruction has historically leaned on per-point tree traversal, the pattern hadn't been applied in my toolchain, and my recognition path ran through a game engine rather than the GPU literature. If the Minecraft framing helps you spot the same shape in your own pipeline, it did its job. The GPU-computing community got there first, independently, and without any game-engine lineage.

The optimization stack

Four optimizations, applied sequentially:

Stage 1: Rust + Rayon (4.1× over scipy). Port the per-voxel solve from Python to Rust with Rayon across 32 cores. Same KD-tree, same kernel, same linear algebra. The gain is the well-known interpreter-overhead elimination and nothing more.

Stage 2: Chunk-based GPU RBF (344× over per-voxel Rust). The architectural change described above: partition the 128³ domain into chunks, run one shared KD-tree lookup per chunk, and batch all of a chunk's per-voxel systems into a single GPU call (NVIDIA RTX 4090). Independent solves drop from about 2 million to about 4,096.

Stage 3: Gram-matrix eigendecomposition (2.2× over full SVD). The matrix-completion solver in the same pipeline (Soft-Impute) computed a full SVD of a 16,384 × 128 matrix on each iteration. For a tall-and-skinny matrix, the singular values come cheaper from the eigendecomposition of the 128 × 128 Gram matrix MᵀM. This is a standard shortcut.

Stage 4: Dual-machine parallelism (~2×). Independent experiments distributed across two machines (an RTX 4090 workstation and an RTX 3060 Ti server) over a Tailscale mesh.

Each speedup is measured against its own baseline, and the stages target different pipeline components, so they don't meaningfully multiply into one grand total. The honest end-to-end number is this: the validation campaign estimated at 57 hours under the naive implementation measured 17.6 minutes after the full stack, a 194× improvement in practice.

What the numbers do and don't say

A few things the review process forced me to be precise about:

The 344× bundles hardware and architecture. Stage 2 changes two things at once: the algorithmic structure (chunked shared lookups with batched solves instead of per-voxel tree traversal) and the execution hardware (an RTX 4090 instead of a 32-core CPU). I did not run a GPU-side spatial-index baseline such as a GPU KD-tree or BVH, so I can't say how much of the 344× comes from the representation change and how much comes from the well-documented gain of moving any spatial query to a GPU. The literature warns explicitly about conflating the two (Lee et al. 2010), and GPU tree traversal alone is typically worth 10 to 30× over CPU. A fair decomposition would likely attribute a large share to the hardware. The before-and-after is real, but the attribution is unmeasured.

Numerical agreement, stated exactly. The chunked GPU implementation was validated against the reference implementation. Cross-machine runs, covering WSL2 versus native Linux, Ada versus Ampere GPU architectures, and different CUDA runtimes, agreed to three decimal places. I'm not claiming bit-identical output; batched GPU linear algebra doesn't work that way.

Scope. All of this is specific to RBF interpolation on regular voxel grids with uniform chunk decomposition. Irregular meshes, adaptive grids, or different kernels may not benefit. GPU memory bounds the chunk size, so smaller cards mean smaller chunks and more kernel launches.

What it actually changed

The speedup itself is not the interesting outcome. The interesting outcome is what became affordable.

Before, one reconstruction was a 4.3-hour commitment, and proper multi-seed validation was the kind of thing you plan a week around, which in practice means it doesn't happen. After, the full four-experiment adversarial campaign of 380 complete reconstructions across two machines finished in 25.8 minutes of wall clock. Multi-seed statistical validation went from too expensive to do properly to something you run while the coffee brews. Every adversarial robustness result in the density study exists because this restructuring made the campaign cheap enough to run at all, and the same RBF workload also underpins the earlier ocean reconstruction study.

The takeaway is an old one that the GPU literature has understood for nearly two decades: how you arrange a computation for the hardware, meaning data locality, batching boundaries, and memory movement, can matter more than what language or library you write it in. The research companion to this post explores a statistical cousin of that idea (representation-dependent recoverability). On the engineering side there is no new principle on offer here, just a reminder that boring-sounding restructuring is sometimes worth two orders of magnitude, and that useful patterns can arrive by unexpected routes. Mine arrived by way of Minecraft.

A note on how this was built

I didn't write the implementation code. Claude Code (Anthropic) implemented the optimization stack from my architectural specifications, Claude managed research state across sessions, and ChatGPT (OpenAI) provided independent adversarial review, including the prior-art review that correctly demoted this from paper to post. The chunk-decomposition strategy, the GPU batching approach, and the dual-machine validation protocol were my calls. The line-by-line code was not.

References

• Anderson, J. A., Lorenz, C. D., & Travesset, A. (2008). General purpose molecular dynamics simulations fully implemented on graphics processing units. Journal of Computational Physics, 227(10).
• Green, S. (2008). Particle Simulation using CUDA. NVIDIA Technical Report (CUDA SDK).
• Hardy, D. J., Stone, J. E., & Schulten, K. (2009). Multilevel summation of electrostatic potentials using graphics processing units. Parallel Computing, 35(3).
• Kirk, D. B., & Hwu, W. W. (2011). Programming Massively Parallel Processors: A Hands-on Approach. Morgan Kaufmann. (The "Input Binning" pattern.)
• Lee, V. W., et al. (2010). Debunking the 100X GPU vs. CPU myth: an evaluation of throughput computing on CPU and GPU. ISCA 2010.
• Rodrigues, C. I., Hardy, D. J., Stone, J. E., Schulten, K., & Hwu, W. W. (2008). GPU acceleration of cutoff pair potentials for molecular modeling applications. ACM Computing Frontiers 2008.
• Poehlein, J. A. (2026). Same Field, Opposite Winners: Representation-Dependent Recoverability in Spatial Reconstruction. mukenshi.com/research/cosmic-density-reconstruction.
• Poehlein, J. A. (2026). Can Netflix Math Map the Ocean? A Negative Result Worth Sharing. mukenshi.com/research/ocean-matrix-completion.