GPTQ: Mathematical Foundations, Optimal Brain Surgeon Inversion, and Second-Order Error Minimization in LLM Quantization

Large language model inference during autoregressive generation is overwhelmingly memory-bandwidth bound. For batch size 1 decoding, each generated token requires streaming every parameter of a model from High Bandwidth Memory (HBM) into GPU SRAM and Tensor Cores. A 70-billion parameter model in 16-bit precision (FP16 or BF16) requires roughly 140 GB of VRAM, exceeding the capacity of a single 80 GB NVIDIA A100 or H100 GPU and demanding multi-GPU tensor parallelism solely to hold the model weights.

Post-training quantization (PTQ) compresses model weights into low-bit integer formats (such as INT4 or INT3) without requiring full retraining. However, naive round-to-nearest (RTN) quantization causes severe perplexity degradation below 8 bits due to complex cross-weight interactions and activation outliers.

GPTQ (Accurate Post-Training Quantization for Generative Pre-trained Transformers), introduced by Elias Frantar, Saleh Ashkboos, Torsten Hoefler, and Dan Alistarh (ICLR 2023), resolved this problem. GPTQ adapts the classical second-order Optimal Brain Surgeon framework into an efficient, hardware-aware algorithm capable of quantizing 175-billion parameter models in approximately 4 GPU hours with negligible loss in accuracy.

This explainer details the mathematical derivation of second-order weight compensation, the algorithmic innovations that make GPTQ computationally tractable, its numerical stabilization techniques, and its serving mechanics.

The Layer-Wise Reconstruction Problem

Post-training quantization operates layer by layer to avoid the intractable computational cost of optimizing all model parameters simultaneously over large corpora.

Let a linear layer have a full-precision weight matrix $\mathbf{W} \in \mathbb{R}^{d_{\text{out}} \times d_{\text{in}}}$ and receive a calibration input matrix $\mathbf{X} \in \mathbb{R}^{d_{\text{in}} \times N}$, where $N$ is the total number of calibration tokens passed through the network. The goal is to find a quantized weight matrix $\widehat{\mathbf{W}} \in \mathcal{Q}^{d_{\text{out}} \times d_{\text{in}}}$ that minimizes the squared error of the output activations:

$$\min_{\widehat{\mathbf{W}}} \|\mathbf{W}\mathbf{X} - \widehat{\mathbf{W}}\mathbf{X}\|_F^2$$

where $\|\cdot\|_F$ denotes the Frobenius norm and $\mathcal{Q}$ represents the discrete quantization grid (for example, signed or unsigned integers scaled by a step size).

Because the Frobenius norm is the sum of squared Euclidean norms of the matrix rows, the layer-wise objective decomposes into $d_{\text{out}}$ independent subproblems:

$$\|\mathbf{W}\mathbf{X} - \widehat{\mathbf{W}}\mathbf{X}\|_F^2 = \sum_{i=1}^{d_{\text{out}}} \|\mathbf{w}_{i,:} \mathbf{X} - \widehat{\mathbf{w}}_{i,:} \mathbf{X}\|_2^2$$

where $\mathbf{w}_{i,:} \in \mathbb{R}^{1 \times d_{\text{in}}}$ is the $i$-th row of $\mathbf{W}$. Each row can therefore be optimized independently using identical input statistics.

       Input Activations X (d_in x N)
                 │
                 ▼
┌─────────────────────────────────────────────────┐
│ Full-Precision Weights W  (d_out x d_in)        │
└───────────────────────┬─────────────────────────┘
                        │
                        ▼ Layer-Wise Error Minimization
┌─────────────────────────────────────────────────┐
│ Quantized Weights W_hat   (d_out x d_in)        │
│ s.t. argmin || W X - W_hat X ||_F^2             │
└─────────────────────────────────────────────────┘

Second-Order Taylor Expansion and the Hessian

Consider a single row $\mathbf{w} \in \mathbb{R}^{d_{\text{in}}}$ (written as a column vector for notation convenience) and an arbitrary perturbation vector $\boldsymbol{\delta} \in \mathbb{R}^{d_{\text{in}}}$ representing the quantization error $\boldsymbol{\delta} = \widehat{\mathbf{w}} - \mathbf{w}$.

The squared output error $E(\mathbf{w} + \boldsymbol{\delta}) = \frac{1}{2} \|(\mathbf{w} + \boldsymbol{\delta})^T \mathbf{X} - \mathbf{w}^T \mathbf{X}\|_2^2$ can be expanded via a Taylor series around the unquantized weights $\mathbf{w}$:

$$E(\mathbf{w} + \boldsymbol{\delta}) = E(\mathbf{w}) + \mathbf{g}^T \boldsymbol{\delta} + \frac{1}{2} \boldsymbol{\delta}^T \mathbf{H} \boldsymbol{\delta}$$

Here, $\mathbf{g} = \nabla_{\mathbf{w}} E(\mathbf{w})$ is the gradient and $\mathbf{H} = \nabla^2_{\mathbf{w}} E(\mathbf{w})$ is the Hessian matrix.

Because $\mathbf{w}$ is the exact unquantized weight vector, the reconstruction error at $\boldsymbol{\delta} = \mathbf{0}$ is exactly zero: $E(\mathbf{w}) = 0$. Consequently, the gradient $\mathbf{g}$ at $\boldsymbol{\delta} = \mathbf{0}$ vanishes ($\mathbf{g} = \mathbf{0}$).

The error simplifies strictly to the quadratic curvature term:

$$E(\boldsymbol{\delta}) = \frac{1}{2} \boldsymbol{\delta}^T \mathbf{H} \boldsymbol{\delta}$$

Computing the Hessian with respect to $\mathbf{w}$ yields:

$$\mathbf{H} = 2 \mathbf{X} \mathbf{X}^T \in \mathbb{R}^{d_{\text{in}} \times d_{\text{in}}}$$

The Hessian is twice the unnormalized sample covariance of the input activations. Crucially, $\mathbf{H}$ depends solely on the layer inputs $\mathbf{X}$ and is completely independent of the weight values $\mathbf{w}$. Therefore, the exact same Hessian matrix $\mathbf{H}$ governs the error curvature for every single row $i \in \{1, \dots, d_{\text{out}}\}$ of the weight matrix $\mathbf{W}$.

Optimal Brain Surgeon (OBS) and Optimal Brain Quantization (OBQ)

The framework of using second-order Taylor expansions to modify neural network weights traces back to Optimal Brain Damage (LeCun et al., 1989) and Optimal Brain Surgeon (Hassibi & Stork, 1993) for pruning, later generalized to quantization in Optimal Brain Quantization (Frantar et al., 2022).

Suppose we select a specific weight index $q \in \{1, \dots, d_{\text{in}}\}$ and round it to its quantized value $\text{quant}(w_q)$. We wish to adjust all remaining unquantized weights (the set of free indices $F$) to minimize the total quadratic error.

The problem is formulated as a constrained quadratic program:

$$\min_{\boldsymbol{\delta}} \frac{1}{2} \boldsymbol{\delta}^T \mathbf{H} \boldsymbol{\delta} \quad \text{subject to} \quad \mathbf{e}_q^T \boldsymbol{\delta} = \text{quant}(w_q) - w_q$$

where $\mathbf{e}_q$ is the $q$-th standard basis vector.

Lagrangian Derivation

We define the Lagrangian function with Lagrange multiplier $\lambda \in \mathbb{R}$:

$$\mathcal{L}(\boldsymbol{\delta}, \lambda) = \frac{1}{2} \boldsymbol{\delta}^T \mathbf{H} \boldsymbol{\delta} + \lambda \left( \mathbf{e}_q^T \boldsymbol{\delta} - (\text{quant}(w_q) - w_q) \right)$$

Taking the gradient with respect to $\boldsymbol{\delta}$ and setting it to zero:

$$\nabla_{\boldsymbol{\delta}} \mathcal{L} = \mathbf{H} \boldsymbol{\delta} + \lambda \mathbf{e}_q = \mathbf{0} \implies \boldsymbol{\delta} = -\lambda \mathbf{H}^{-1} \mathbf{e}_q$$

Multiplying both sides by $\mathbf{e}_q^T$ to apply the constraint:

$$\mathbf{e}_q^T \boldsymbol{\delta} = -\lambda \mathbf{e}_q^T \mathbf{H}^{-1} \mathbf{e}_q = -\lambda [\mathbf{H}^{-1}]_{qq} = \text{quant}(w_q) - w_q$$

Solving for the Lagrange multiplier $\lambda$:

$$\lambda = -\frac{\text{quant}(w_q) - w_q}{[\mathbf{H}^{-1}]_{qq}}$$

Substituting $\lambda$ back into the expression for $\boldsymbol{\delta}$:

$$\boldsymbol{\delta} = -\frac{w_q - \text{quant}(w_q)}{[\mathbf{H}^{-1}]_{qq}} \cdot (\mathbf{H}^{-1})_{:, q}$$

where $(\mathbf{H}^{-1})_{:, q}$ is the $q$-th column of the inverse Hessian matrix.

Minimal Error Penalty

Substituting the optimal compensation vector $\boldsymbol{\delta}$ back into the objective function $E(\boldsymbol{\delta}) = \frac{1}{2} \boldsymbol{\delta}^T \mathbf{H} \boldsymbol{\delta}$:

$$\Delta E_q = \frac{1}{2} \left( -\lambda \mathbf{H}^{-1} \mathbf{e}_q \right)^T \mathbf{H} \left( -\lambda \mathbf{H}^{-1} \mathbf{e}_q \right) = \frac{1}{2} \lambda^2 \mathbf{e}_q^T \mathbf{H}^{-1} \mathbf{e}_q = \frac{1}{2} \lambda^2 [\mathbf{H}^{-1}]_{qq}$$

Substituting $\lambda$:

$$\Delta E_q = \frac{(w_q - \text{quant}(w_q))^2}{2 [\mathbf{H}^{-1}]_{qq}}$$

The quantity $\Delta E_q$ defines the second-order saliency of weight $w_q$. It measures the exact increase in reconstruction error incurred when quantizing $w_q$ and optimally compensating all remaining unquantized weights.

Recursive Inverse Hessian Updating

After quantizing weight index $q$, index $q$ is removed from the set of active weights $F$. The inverse Hessian of the remaining submatrix $\mathbf{H}_{F \setminus \{q\}}$ must be updated. By block matrix inversion and Gaussian elimination, the update rule is:

$$\mathbf{H}_{F \setminus \{q\}}^{-1} = (\mathbf{H}_F^{-1})_{F \setminus \{q\}, F \setminus \{q\}} - \frac{1}{[\mathbf{H}_F^{-1}]_{qq}} (\mathbf{H}_F^{-1})_{F \setminus \{q\}, q} (\mathbf{H}_F^{-1})_{q, F \setminus \{q\}}$$

This update eliminates row and column $q$ from $\mathbf{H}_F^{-1}$ in $O(d^2)$ operations via a symmetric rank-1 update.

The Computational Bottleneck of OBQ

While OBQ produces highly accurate quantized networks, it cannot scale directly to large language models.

OBQ Computational Scaling Bottleneck:
  • Greedy selection: Evaluates argmin Delta E_q across all remaining parameters
  • Row-by-row execution: Quantizes each of d_out rows separately
  • Scalar rank-1 updates: Operates at BLAS Level 1/2 (memory bandwidth bound)
  • Complexity: O(d_out * d_in^3) -> Weeks of compute for 175B parameter models

For a standard transformer layer with $d_{\text{in}} = d_{\text{out}} = 12,288$, quantizing a single weight matrix via OBQ requires:

$$\mathcal{O}(d_{\text{out}} \cdot d_{\text{in}}^3) \approx 12,288 \times (12,288)^3 \approx 2.28 \times 10^{16} \text{ FLOPs}$$

Because each scalar step involves small vector operations (BLAS Level 1 and 2), memory access latency dominates GPU execution, utilizing less than 5% of peak Tensor Core throughput. Quantizing an entire 175B model would take multiple weeks.

The GPTQ Algorithmic Innovations

GPTQ overcomes the OBQ computational bottleneck through three distinct algorithmic advancements: arbitrary order quantization, lazy batch updates, and Cholesky reformulation.

1. Arbitrary Order Quantization

The standard OBQ algorithm greedily quantizes the weight with the smallest error $\Delta E_q$ at each step. This greedy choice couples the quantization order to the individual weight values of each row, forcing every row to follow a different quantization trajectory.

Frantar et al. made a critical empirical discovery: quantizing weights in an arbitrary, fixed column order (e.g. simply left-to-right from column $1$ to column $d_{\text{in}}$) achieves virtually identical perplexity on large models compared to greedy selection.

Because the quantization order is fixed across all rows:

• The sequence of inverse Hessian matrices $\mathbf{H}_F^{-1}$ and elimination steps is identical for all $d_{\text{out}}$ rows.
• The inverse Hessian updates need to be computed exactly once per layer, rather than $d_{\text{out}}$ times.
• All $d_{\text{out}}$ rows of $\mathbf{W}$ can be updated simultaneously in batched vector and matrix operations.

2. Lazy Batch Updates (Block Quantization)

Updating the remaining full-precision columns across the entire matrix $\mathbf{W}$ after quantizing every single column $j$ requires repeated memory round-trips. Each column update performs a rank-1 outer product $\mathbf{W}_{:, F} \leftarrow \mathbf{W}_{:, F} - \boldsymbol{\delta}_{:, j} (\mathbf{H}^{-1})_{j, F}$, which is memory-bandwidth bound.

To maximize arithmetic intensity, GPTQ groups columns into blocks of size $B$ (typically $B = 128$). The quantization proceeds in two stages:

1. Intra-Block Updates (Fast Local Updates): Inside a block of $B$ columns, weights are sequentially quantized. After column $j \in B$ is rounded, only the remaining columns within the block $B$ are immediately updated using fast cache-resident operations.
2. Inter-Block Global Updates (BLAS-3 GEMM): The accumulated quantization errors across the entire block, $\mathbf{E}_B = \mathbf{W}_{:, B} - \widehat{\mathbf{W}}_{:, B} \in \mathbb{R}^{d_{\text{out}} \times B}$, are gathered. The remaining unquantized columns $F$ (all columns to the right of block $B$) are updated in a single large Matrix-Matrix Multiplication (GEMM):

$$\mathbf{W}_{:, F} \leftarrow \mathbf{W}_{:, F} - \mathbf{E}_B \cdot \left( [\mathbf{H}^{-1}]_{B, B} \right)^{-1} [\mathbf{H}^{-1}]_{B, F}$$

By offloading the vast majority of FLOPs to high-throughput GEMM operations (BLAS Level 3), GPTQ achieves near-peak GPU hardware utilization.

3. Cholesky Decomposition and Pre-computation

Instead of computing sequential rank-1 updates to $\mathbf{H}^{-1}$ at runtime, GPTQ uses the Cholesky decomposition of the inverse Hessian.

Because $\mathbf{H} = 2\mathbf{X}\mathbf{X}^T$ is symmetric positive semi-definite, its inverse $\mathbf{H}^{-1}$ is also symmetric positive definite. Computing the Cholesky factorization:

$$\mathbf{H}^{-1} = \mathbf{L} \mathbf{L}^T$$

where $\mathbf{L}$ is a lower triangular matrix.

The entries of $\mathbf{L}$ directly encode the unquantized weight compensation weights. Specifically, the required update vector for column $j$ onto remaining columns $k > j$ is given directly by:

$$\frac{(\mathbf{H}_F^{-1})_{k, j}}{[\mathbf{H}_F^{-1}]_{jj}} = \frac{L_{k, j}}{L_{j, j}}$$

and the diagonal element is:

$$[\mathbf{H}_F^{-1}]_{jj} = L_{j, j}^2$$

This formulation eliminates numerical instability from repeated subtractions in rank-1 updates and allows all compensation coefficients to be pre-calculated in a single initial Cholesky step.

Hessian Damping and Numerical Regularization

In deep transformer networks, activations across different hidden dimensions are often highly correlated or contain zero-variance channels. Consequently, the raw activation covariance matrix $\mathbf{X}\mathbf{X}^T$ can be ill-conditioned or numerically singular.

Inversion of an ill-conditioned Hessian leads to numerical explosion in $[\mathbf{H}^{-1}]_{jj}$, causing massive, destructive weight compensation updates.

GPTQ introduces a symmetric diagonal damping term $\lambda \mathbf{I}$:

$$\mathbf{H}_{\text{damped}} = 2 \mathbf{X} \mathbf{X}^T + \lambda \mathbf{I}$$

$$\lambda = \alpha \cdot \frac{1}{d_{\text{in}}} \operatorname{Tr}(\mathbf{H}) = \alpha \cdot \operatorname{mean}(\operatorname{diag}(\mathbf{H}))$$

The damping factor $\alpha$ is typically set to $0.01$ (a 1% ridge penalty).

This regularization provides two mathematical guarantees:

1. It shifts all eigenvalues strictly into the positive domain: $\lambda_{\min}(\mathbf{H}_{\text{damped}}) \ge \lambda > 0$, ensuring the matrix is strictly invertible and well-conditioned for Cholesky decomposition.
2. It bounds the maximum compensation step size, preventing outlier activations from destabilizing neighboring weights.

Activation Ordering (Act-Order / Desc-Act)

Although arbitrary order quantization works well on average, certain transformer layers exhibit extreme sensitivity where specific input channels carry disproportionate variance.

To improve quantization accuracy on sensitive layers, GPTQ supports Activation Ordering (also called act-order or desc_act):

1. Compute the diagonal elements of the Hessian matrix $H_{jj} = 2 \sum_{k=1}^N X_{jk}^2$, which measure the total activation energy of each input channel $j$.
2. Permute the columns of the weight matrix $\mathbf{W}$ and rows/columns of $\mathbf{H}$ such that channels with the largest diagonal values are quantized first:

$$\pi = \operatorname{argsort}\left( \operatorname{diag}(\mathbf{H}), \text{descending} \right)$$

1. Apply GPTQ quantization on the permuted matrix $\mathbf{W}_{:, \pi}$.
2. Invert the permutation on the quantized weights $\widehat{\mathbf{W}}$ to restore the original tensor layout.

By quantizing the highest-energy channels first, the algorithm leaves the maximum number of remaining unquantized parameters ($F$) available to absorb and compensate for the errors in those critical channels.

Act-Order Workflow:
  1. Compute diagonal activation energy: diag(H) = 2 * sum(X^2, dim=tokens)
  2. Sort indices by energy: pi = argsort(diag(H), descending=True)
  3. Permute weights & Hessian: W_sorted = W[:, pi], H_sorted = H[pi, :][:, pi]
  4. Run Block-Wise GPTQ on W_sorted
  5. Unpermute quantized weights: W_quant = W_quant[:, inv_pi]

Quantization Formats and Group Size

Quantization maps continuous weights $w$ to discrete integers $q \in \{0, \dots, 2^b - 1\}$ using a scale $s$ and zero-point $z$:

$$\widehat{w} = s \cdot (q - z)$$

$$q = \operatorname{clamp}\left( \left\lfloor \frac{w}{s} \right\rceil + z, 0, 2^b - 1 \right)$$

For symmetric quantization, $z = 0$ (or $2^{b-1}$), and the grid is centered at zero:

$$s = \frac{\max(|w|)}{2^{b-1} - 1}$$

Per-Channel vs. Group-Wise Quantization

• Per-Channel (Per-Column / Per-Row): A single scale $s$ and zero-point $z$ are assigned to an entire output channel ($1 \times d_{\text{in}}$).
• Group-Wise Quantization: Each row is split into independent sub-vectors of size $g$ (typically $g = 128$ or $g = 64$), each with its own scale $s_k$ and zero-point $z_k$.

Group-wise quantization dramatically reduces the dynamic range within any single quantization bucket, virtually eliminating the impact of localized weight outliers at the cost of modest metadata storage overhead:

$$\text{Bit-rate per parameter} = b + \frac{16 + 16}{g} \text{ bits}$$

For $b = 4$ and $g = 128$, the metadata overhead is $\frac{32}{128} = 0.25$ bits per weight, yielding an effective bit-rate of 4.25 bits per parameter.

GPTQ Algorithm Execution Walkthrough

The following PyTorch-style algorithm outlines the complete GPTQ process with lazy batch updates:

import torch

def gptq_quantize_layer(
    W: torch.Tensor,          # Shape: [d_out, d_in] (FP16/FP32)
    H: torch.Tensor,          # Shape: [d_in, d_in] (2 * X @ X.T)
    bits: int = 4,
    block_size: int = 128,
    percdamp: float = 0.01,
    group_size: int = 128
) -> torch.Tensor:
    d_out, d_in = W.shape
    W_quant = W.clone().float()
    
    # 1. Hessian Damping Regularization
    dead = torch.diag(H) == 0
    H[dead, dead] = 1.0
    W_quant[:, dead] = 0.0
    
    diag_mean = torch.mean(torch.diag(H))
    damp = percdamp * diag_mean
    H += damp * torch.eye(d_in, device=H.device)
    
    # 2. Inversion via Cholesky Decomposition
    H_inv = torch.linalg.cholesky_inverse(torch.linalg.cholesky(H))
    H_inv_chol = torch.linalg.cholesky(H_inv, upper=True)
    
    # 3. Process columns in blocks
    for block_start in range(0, d_in, block_size):
        block_end = min(block_start + block_size, d_in)
        count = block_end - block_start
        
        W_block = W_quant[:, block_start:block_end].clone()
        H_inv_block = H_inv[block_start:block_end, block_start:block_end]
        H_inv_chol_block = H_inv_chol[block_start:block_end, block_start:block_end]
        
        Err_block = torch.zeros_like(W_block)
        
        # Intra-block quantization loop
        for j in range(count):
            col_idx = block_start + j
            w = W_block[:, j]
            d = H_inv_chol_block[j, j]
            
            # Scalar quantization (per-channel or group-wise)
            q = quantize_weights(w, bits=bits, group_size=group_size)
            Err_block[:, j] = (w - q) / d
            
            # Update remaining columns within the local block
            W_block[:, j:] -= Err_block[:, j:j+1] @ H_inv_chol_block[j:j+1, j:]
        
        W_quant[:, block_start:block_end] = quantize_weights(
            W_quant[:, block_start:block_end], bits=bits, group_size=group_size
        )
        
        # Inter-block global update via BLAS-3 GEMM
        if block_end < d_in:
            H_inv_remaining = H_inv[block_start:block_end, block_end:]
            W_quant[:, block_end:] -= Err_block @ H_inv_remaining

    return W_quant

Comparison: GPTQ vs. Alternative Quantization Schemes

Method Comparison

• RTN (Round-to-Nearest): Weight-only (W4A16). Heuristic rounding without optimization. Calibration runs in seconds with minimal memory overhead, but lacks outlier handling and degrades catastrophically below 8-bit precision.
• GPTQ: Weight-only (W4A16 / W3A16 / W2A16). Second-order inverse Hessian error minimization. Calibrates a 70B parameter model in 15 to 30 minutes. Delivers very high hardware throughput via Marlin and ExLlama kernels while compensating errors across remaining weights.
• AWQ: Weight-only (W4A16). First-order activation-aware per-channel scaling. Calibrates a 70B model in 15 to 30 minutes. Protects salient weights by scaling input channels and achieves high throughput via W4A16 GEMM kernels.
• SmoothQuant: Weight and activation quantization (W8A8). Exact mathematical transformation ($Y = (X \cdot S^{-1}) \cdot (S \cdot W)$). Calibrates in minutes, migrating activation outlier difficulty into weights to enable standard INT8 Tensor Core GEMM.
• QuIP# / AQLM: Vector quantization (2-bit). Employs randomized Hadamard transforms, incoherence processing, and vector codebooks. Calibration takes 2 to 6 hours, trading codebook lookup complexity for high fidelity at 2-bit precision.

GPTQ vs. AWQ

While GPTQ optimizes weight rounding through second-order inverse Hessian error compensation, AWQ (Activation-aware Weight Quantization; Lin et al., 2023) observes that not all weights are equally important. AWQ protects the top 1% of weights corresponding to high-magnitude activation channels by applying per-channel scale factors $S$ to activations and inverse scaling $S^{-1}$ to weights before applying standard rounding.

In practice, both methods achieve comparable perplexity at 4-bit precision. GPTQ provides superior flexibility for extreme low-bit regimes (such as 3-bit and 2-bit quantization) because its Hessian compensation actively redistributes quantization noise across all parameters.

Inference Execution and Modern Kernel Architectures

During inference, weights stored in 4-bit packed integer formats (e.g., eight 4-bit weights packed into a single int32 word) must be unpacked and dequantized to FP16/BF16 prior to matrix multiplication.

Dequantization Kernels

1. ExLlama (v1 / v2): Implements fused dequantize-and-GEMV CUDA kernels tailored for low-batch autoregressive generation (batch sizes 1-4). ExLlama streams packed INT4 weights directly into GPU registers, performs fast bit-shifting and scale-multiplication in register space, and accumulates FP16 dot products directly on CUDA cores.
2. Marlin (Mixed Auto-Regressive Linear): A highly optimized FP16xINT4 GEMM kernel designed by Elias Frantar et al. (2024). Marlin restructures memory layouts into 2D tiles to exploit both memory bandwidth and Tensor Core asynchronous copy instructions (cp.async), delivering near-ideal linear speedups across both low-batch and high-batch inference regimes.
3. FlashInfer: Provides composable GPU kernel templates that support group-quantized GPTQ formats with paged KV caching and unified batching schedules.

FP16xINT4 Fused Inference Pipeline:
  1. Load packed INT4 weights from HBM (4x memory bandwidth reduction)
  2. Stream FP16 activations from HBM/SRAM
  3. Fused register-level bit-shift and dequantization (W_fp16 = (W_int4 - zero) * scale)
  4. Tensor Core GEMM accumulation in FP16 / FP32
  5. Output FP16 activation tensor

Summary

GPTQ transformed post-training quantization from an impractical theoretical framework into an industry-standard compression technique for large language models.

By framing layer-wise quantization through a second-order Taylor expansion and proving that arbitrary column ordering allows uniform inverse Hessian updates across all rows, GPTQ scaled the classical Optimal Brain Surgeon formulation to models with hundreds of billions of parameters. Combined with lazy block updates, Cholesky factorization, diagonal damping, and modern fused dequantization kernels, GPTQ enables 4-bit and 3-bit LLM serving with near-lossless output fidelity and substantial memory bandwidth reductions.

Sources

• Frantar, E., Ashkboos, S., Hoefler, T., & Alistarh, D. (2022). GPTQ: Accurate Post-Training Quantization for Generative Pre-trained Transformers. arXiv:2210.17323. https://arxiv.org/abs/2210.17323
• Hassibi, B., & Stork, D. G. (1993). Second order derivatives for network pruning: Optimal Brain Surgeon. Advances in Neural Information Processing Systems (NeurIPS 1992). https://proceedings.neurips.cc/paper/1992/hash/303ed4c69846ab36c2904d3ba8573050-Abstract.html
• LeCun, Y., Denker, J. S., & Solla, S. A. (1989). Optimal Brain Damage. Advances in Neural Information Processing Systems (NeurIPS 1989). https://proceedings.neurips.cc/paper/1989/hash/6c9882bbac1c7093bd25041881277658-Abstract.html
• Frantar, E., & Alistarh, D. (2022). Optimal Brain Compression: A Framework for Accurate Post-Training Quantization and Pruning. Advances in Neural Information Processing Systems (NeurIPS 2022). https://arxiv.org/abs/2208.11580
• Lin, J., Tang, J., Tang, H., Yang, S., Dang, X., & Han, S. (2023). AWQ: Activation-aware Weight Quantization for LLM Compression and Acceleration. arXiv:2306.00978. https://arxiv.org/abs/2306.00978
• Xiao, G., Lin, J., Mickelson, M., Han, S., & Tang, S. (2023). SmoothQuant: Accurate and Efficient Post-Training Quantization for Large Language Models. ICML 2023. arXiv:2211.10438. https://arxiv.org/abs/2211.10438
• Frantar, E., & Alistarh, D. (2024). Marlin: A Fast FP16xINT4 Matrix Multiplication Kernel for Efficient LLM Inference. arXiv:2408.11743. https://arxiv.org/abs/2408.11743