flash-sampling-efficient-decoding

name	flash-sampling-efficient-decoding
title	FlashSampling: Fast and Memory-Efficient Exact Sampling
version	0.0.2
engine	skillxiv-v0.0.2-claude-opus-4.6
license	MIT
url	https://arxiv.org/abs/2603.15854
keywords	["Sampling Efficiency","Exact Sampling","Gumbel Trick","Bandwidth Optimization","Token Generation"]
description	Fuse categorical sampling directly into LM-head matrix multiplication to eliminate logits materialization. Use Gumbel noise during computation and hierarchical reduction to achieve 19% token-level latency reduction.

FlashSampling: Fusing Sampling into Matrix Multiplication

Sampling from the output logits is a hidden performance bottleneck in language model inference. Current approaches materialize large logits tensors in high-bandwidth memory (HBM), consuming significant memory bandwidth and requiring multiple GPU kernels after LM-head computation. FlashSampling solves this by fusing sampling directly into the LM-head matrix multiplication: logits are computed tile-by-tile on-chip, Gumbel noise is applied during computation, and only the maximum per-tile is tracked—enabling exact categorical sampling without materializing full logits. This achieves up to 19% reduction in per-token latency across modern GPUs.

The approach remains mathematically exact because argmax decomposes over partitions (max of maximums equals global maximum).

Core Concept

FlashSampling exploits the mathematical property that sampling from a categorical distribution can be reformulated as:

Standard Sampling:

logits = LM_head(hidden)                    # [vocab_size]
probs = softmax(logits)
token = categorical_sample(probs)

FlashSampling:

# Compute logits tile-by-tile, maintaining only running maximum
max_logit = -inf
max_idx = -1

for tile_idx in range(num_tiles):
    tile_logits = LM_head_partial(hidden, tile_idx)
    tile_logits += Gumbel_noise
    tile_max_idx = argmax(tile_logits)
    tile_max = tile_logits[tile_max_idx]

    if tile_max > max_logit:
        max_logit = tile_max
        max_idx = tile_idx * tile_size + tile_max_idx

# Token with highest Gumbel-perturbed logit is sampled token
token = max_idx

Because argmax commutes with Gumbel perturbation, this is mathematically equivalent to exact sampling.

Architecture Overview

• LM-Head Tiling — Decompose matrix multiplication into independent tiles
• On-Chip Computation — Keep tiles in fast cache, avoid HBM spills
• Gumbel Noise Injection — Apply noise during tile computation
• Running Maximization — Track (value, index) of best tile element
• Grouped Variant — Hierarchical reduction for tensor parallelism
• Kernel Fusion — Single GPU kernel replaces compute + memory stages
• Batch Processing — Maintain efficiency across batch dimensions

Implementation Steps

Start by implementing the basic tile-based sampling logic.

import torch
import torch.nn.functional  F
 numpy  np

 :
    

     ():
        .vocab_size = vocab_size
        .hidden_dim = hidden_dim
        .tile_size = tile_size
        .num_tiles = (vocab_size + tile_size - ) // tile_size


     ():
        
        
        u = torch.rand_like(logits)
        gumbel_noise = -torch.log(-torch.log(u + ) + )

        
        perturbed = (logits + gumbel_noise) / temperature

        
         torch.argmax(perturbed, dim=-)

     () -> torch.Tensor:
        
        batch_size = hidden.size()
        best_indices = torch.zeros(batch_size, dtype=torch.long)
        best_scores = torch.full((batch_size,), ())

        
         tile_idx  (.num_tiles):
            start_vocab = tile_idx * .tile_size
            end_vocab = (start_vocab + .tile_size, .vocab_size)

            
            tile_weights = lm_head_weight[start_vocab:end_vocab, :]
            

            
            tile_logits = torch.matmul(hidden, tile_weights.t())
            

            
            u = torch.rand_like(tile_logits)
            gumbel = -torch.log(-torch.log(u + ) + )

            
            tile_scores = (tile_logits + gumbel) / temperature
            

            
            tile_best_scores, tile_best_local_idx = torch.(tile_scores, dim=)

            
            improved = tile_best_scores > best_scores
            best_scores[improved] = tile_best_scores[improved]
            best_indices[improved] = (
                start_vocab + tile_best_local_idx[improved]
            )

         best_indices

     () -> torch.Tensor:
        
        
        group_vocab_size = .vocab_size // group_size

        best_indices = torch.zeros(hidden.size(), dtype=torch.long)
        best_scores = torch.full((hidden.size(),), ())

         group_idx  (group_size):
            
            group_start = group_idx * group_vocab_size
            group_end = (group_idx + ) * group_vocab_size
            group_weights = lm_head_weight[group_start:group_end, :]

            
             tile_in_group  ((group_vocab_size + .tile_size - )
                                      // .tile_size):
                tile_start = group_start + tile_in_group * .tile_size
                tile_end = (tile_start + .tile_size, group_end)

                tile_weights = lm_head_weight[tile_start:tile_end, :]
                tile_logits = torch.matmul(hidden, tile_weights.t())

                
                u = torch.rand_like(tile_logits)
                gumbel = -torch.log(-torch.log(u + ) + )
                tile_scores = (tile_logits + gumbel) / temperature

                
                tile_best_scores, tile_best_idx = torch.(tile_scores, dim=)
                improved = tile_best_scores > best_scores
                best_scores[improved] = tile_best_scores[improved]
                best_indices[improved] = tile_start + tile_best_idx[improved]

         best_indices