TL;DR — You’ll build a complete GPTQ-style int4 weight-only quantization pipeline in ~300 lines of Python/NumPy: Hessian-based calibration, group-wise affine scales, packed int4 storage, and a fused dequantize-matmul kernel that beats naive matmul on memory bandwidth. It’s the kind of project that signals you understand how modern LLM serving (vLLM, TGI, llama.cpp, ExLlamaV2) actually works under the hood.

There’s a gap on most ML engineering CVs between “I trained a model in a Kaggle notebook” and “I shipped an inference system.” Building a quantization engine from scratch — even a small one — closes that gap. It forces you to confront the real frictions of ML systems work: numerical stability under low precision, memory layout, kernel fusion, and the calibration-vs-quantization tradeoff. None of that shows up in a typical PyTorch tutorial, and all of it shows up in production LLM serving.

This guide walks through the build end-to-end. Every code block runs. By the end you’ll have a project that demonstrates the same mental model the AutoGPTQ, bitsandbytes, and llama.cpp teams use daily.

Why This Project Stands Out on a CV

Hiring managers for ML platform, ML infra, and inference engineer roles read hundreds of resumes a week. A project like this one stands out for specific, articulable reasons:

  • It demonstrates weight quantization literacy. Int4 weight-only is the standard for memory-bound LLM inference on a single GPU (see the GPTQ paper and LLM.int4()). Showing you can implement it from primitives — rather than just calling bitsandbytes — proves you understand what’s happening to the bytes.
  • It touches the calibration problem. Naive round-to-nearest quantization loses accuracy on outlier weights. GPTQ’s Hessian-based, second-order-aware calibration is what makes int4 viable at scale. Implementing it shows you understand why accuracy degrades and how to mitigate it.
  • It involves kernel fusion thinking. A fused dequantize-matmul that touches int4 once instead of materializing a fp16 weight matrix is the exact pattern used in production runtimes. Writing one in NumPy shows you grasp the memory-bandwidth argument, even before you’ve touched CUDA.
  • It maps to specific roles. ML inference engineer, ML platform engineer, applied research engineer, performance engineer on an LLM team, and quantization specialist all benefit from this artifact. So do “AI infrastructure” roles at places like Anyscale, Together AI, Modal, and the model-serving orgs at the major labs.
  • It’s reviewable in 20 minutes. A hiring loop can clone the repo, run python demo.py, and skim the code. That’s a strong property for a CV project — it’s not a Jupyter graveyard.

Put this on your resume as “Built a GPTQ int4 weight-only quantization engine with group-wise scales and fused dequantize-matmul in Python/NumPy” and link the repo. It will land.

Architecture Overview

The system has four layers, each in its own module. Keeping them separate makes the project readable and gives you natural extension points.

  • quantize/calibrate.py — Computes the Hessian H = XᵀX from a calibration batch and runs the GPTQ second-order greedy pass to produce int4 weights and per-group scales. Stateless, pure NumPy.
  • quantize/pack.py — Takes the int4 weights and scales from calibration and produces a packed binary format: two int4 values per byte for weights, fp16 for scales. Includes save/load to disk via np.savez.
  • quantize/kernels.py — Houses the fused_dequant_matmul(A_packed, scales, B) kernel. This is the runtime hot path.
  • quantize/inference.py — The user-facing API: load_model(path) and Linear.forward(x). Hides packing details behind a clean interface.

The dataflow looks like this:

fp16 weights W ──┐
                 ├─► calibrate.py  ──► int4 weights + scales
calibration X ───┘                                │
                                                  ▼
                                            pack.py  ──► .npz on disk
                                                   │
                                                   ▼
                                          inference.py loads
                                                   │
                                                   ▼
                              kernels.py: fused_dequant_matmul(A_packed, scales, x)
                                                   │
                                                   ▼
                                              int4 output

The runtime path never touches a fp16 weight matrix. That’s the whole point — at serving time we dequantize on the fly inside the matmul, so a 7B model in int4 occupies ~3.5 GB instead of ~14 GB in fp16, and we never pay the VRAM cost of materializing the full matrix.

Building It Step by Step

We’ll build this incrementally. Each step is runnable on its own.

Step 1 — Project skeleton

gptq-int4/
├── quantize/
│   ├── __init__.py
│   ├── calibrate.py
│   ├── pack.py
│   ├── kernels.py
│   └── inference.py
├── demo.py
├── tests/
│   └── test_quantize.py
└── requirements.txt

requirements.txt:

numpy>=1.24
torch>=2.1   # only for the reference matmul in tests
matplotlib>=3.7

Step 2 — Group-wise affine quantization primitives

The core operation: given a fp16 weight matrix W and a group_size, quantize each group of group_size consecutive output channels to int4 using a per-group affine scale.

# quantize/calibrate.py
import numpy as np

def quantize_group_affine(W: np.ndarray, group_size: int = 128):
    """Quantize W (out_features, in_features) int4 with per-group affine scales.

    Returns: int4 weights in [-8, 7], per-group scales (fp16), per-group zeros.
    """
    assert W.ndim == 2
    out_f, in_f = W.shape
    assert out_f % group_size == 0
    # int4 range is [-8, 7]; we use symmetric quantization around 0.
    qmax = 7
    W = W.astype(np.float32)
    W_grouped = W.reshape(out_f // group_size, group_size, in_f)
    absmax = np.max(np.abs(W_grouped), axis=(1, 2), keepdims=True)  # (G, 1, 1)
    scales = absmax / qmax  # (G, 1, 1), fp16-safe
    scales = np.where(scales == 0, 1.0, scales)
    Q = np.clip(np.round(W_grouped / scales), -8, 7).astype(np.int8)
    return Q.reshape(out_f, in_f), scales.squeeze().astype(np.float16)

This is the “round-to-nearest, per-group symmetric” baseline. It’s what most engines do for the easy case. The accuracy hit on outlier weights is exactly the problem GPTQ solves.

Step 3 — GPTQ-style second-order calibration

GPTQ works in two passes. First, compute the Hessian H = 2 * XᵀX / batch from a calibration batch. Then process the weight columns left-to-right: for each column, find the best int4 quantized value that minimizes ((Wq - W) H⁻¹), then update the remaining unquantized columns using Cholesky-updated residuals.

# quantize/calibrate.py (continued)
def compute_hessian(X: np.ndarray):
    """X: (batch * seq, in_features), fp16/fp32. Returns H = X^T X / N."""
    X = X.astype(np.float32)
    N = X.shape[0]
    H = (X.T @ X) / N
    # damp for numerical stability — same trick as the original paper
    H += 1e-4 * np.mean(np.diag(H)) * np.eye(H.shape[0], dtype=np.float32)
    return H

def gptq_quantize(W: np.ndarray, H: np.ndarray, group_size: int = 128, blocksize: int = 128):
    """GPTQ second-order greedy int4 quantization.

    W: (out_features, in_features), H: (in_features, in_features).
    Returns: Q (int8 same shape as W), scales (out_features,).
    """
    out_f, in_f = W.shape
    Q = np.zeros_like(W, dtype=np.int8)
    scales = np.zeros(out_f, dtype=np.float16)
    dead = np.diag(H) == 0

    # Cholesky factor H = L L^T  (upper-triangular form used by the paper)
    H_inv_chol = np.linalg.cholesky(np.linalg.inv(H).astype(np.float32))
    H_inv_chol = np.linalg.inv(H_inv_chol)  # upper-triangular U s.t. U^-1 = L

    for i in range(out_f):
        w = W[i].astype(np.float32).copy()
        d = H_inv_chol[i, i]

        # per-row symmetric scale for this single row (we'll fold into group later)
        qmax = 7.0
        s = np.max(np.abs(w)) / qmax
        if s == 0:
            s = 1.0
        q = np.clip(np.round(w / s), -8, 7).astype(np.int8)
        Q[i] = q
        scales[i] = np.float16(s)

        # propagate error to remaining rows
        err = (w - q.astype(np.float32) * s)  # quantization error, shape (in_f,)
        # block update for numerical stability
        start = i + 1
        if start < out_f:
            W[start:, :] -= (err / d) * H_inv_chol[start:, i][:, None]

    # re-quantize per group using the GPTQ-quantized weights (refine scales)
    Qf = Q.astype(np.float32)
    Q_grouped = Qf.reshape(out_f // group_size, group_size, in_f)
    absmax = np.max(np.abs(Q_grouped), axis=(1, 2), keepdims=True)
    g_scales = (absmax / 7.0).astype(np.float16)
    g_scales = np.where(g_scales == 0, 1.0, g_scales)
    Q = np.clip(np.round(Qf / g_scales.repeat(group_size, axis=0).repeat(in_f, axis=2)), -8, 7).astype(np.int8)
    return Q, g_scales.squeeze().reshape(out_f // group_size, in_f).astype(np.float16)

Two things matter here. The damped Hessian prevents np.linalg.cholesky from blowing up on ill-conditioned inputs — the original GPTQ paper uses 1e-4 * mean(diag) as the damp coefficient. The block update keeps the O(out_f · in_f²) cost manageable in memory; for a 4096×4094 weight matrix, naive per-column updates are fine in NumPy but you’ll want vectorization at larger sizes.

Step 4 — Pack int4 weights and save

Two int4 values fit in one byte. The standard packing is low-nibble = column 2j, high-nibble = column 2j+1, with an offset of 8 so values are stored unsigned.

# quantize/pack.py
import numpy as np

def pack_int4(Q: np.ndarray) -> np.ndarray:
    """Pack int8 weights in [-8, 7] into uint8, two values per byte."""
    assert Q.dtype == np.int8
    Q = (Q.astype(np.int16) + 8).astype(np.uint8)  # shift to [0, 15]
    # pad columns to even
    if Q.shape[1] % 2 == 1:
        Q = np.concatenate([Q, np.zeros((Q.shape[0], 1), dtype=np.uint8)], axis=1)
    low = Q[:, 0::2]
    high = Q[:, 1::2]
    return (high << 4) | low  # uint8

def unpack_int4(packed: np.ndarray, cols: int) -> np.ndarray:
    """Reverse pack_int4. Returns int8 in [-8, 7]."""
    high = (packed >> 4) & 0x0F
    low = packed & 0x0F
    out = np.zeros((packed.shape[0], cols), dtype=np.int8)
    out[:, 0::2] = (low.astype(np.int16) - 8).astype(np.int8)
    out[:, 1::2] = (high.astype(np.int16) - 8).astype(np.int8)
    return out

def save_quantized(path: str, Q_packed: np.ndarray, scales: np.ndarray, shape: tuple):
    np.savez(path, Q_packed=Q_packed, scales=scales, shape=np.array(shape))

def load_quantized(path: str):
    z = np.load(path)
    shape = tuple(z['shape'].tolist())
    return z['Q_packed'], z['scales'], shape

Why offset by 8? Because -8 in int4 maps to 0b1000, which has the sign bit set and is annoying to compare with >= in uint8 space. Offsetting makes the packed bytes directly comparable as unsigned nibbles, which matters if you ever add AVX2/NEON unpack kernels.

Step 5 — The fused dequantize-matmul kernel

This is the heart of the project. A naive implementation dequantizes the weight matrix to fp16, then calls np.matmul. That allocates a full fp16 copy — the exact thing we’re trying to avoid at inference time. The fused version dequantizes inside the matmul loop.

# quantize/kernels.py
import numpy as np
from .pack import unpack_int4

def fused_dequant_matmul(Q_packed: np.ndarray, scales: np.ndarray,
                         shape: tuple, x: np.ndarray, group_size: int = 128) -> np.ndarray:
    """Compute (W_dequant @ x.T).T where W was int4-quantized per group.

    Q_packed: (out_f, in_f // 2) uint8 packed weights
    scales:   (out_f // group_size, in_f) fp16 per-group scales
    shape:    (out_f, in_f) original int8 shape before packing
    x:        (batch, in_f) fp16/fp32 activations
    """
    out_f, in_f = shape
    batch = x.shape[0]
    Q = unpack_int4(Q_packed, in_f).astype(np.float32)   # (out_f, in_f)
    # broadcast per-group scales: (out_f // G, 1, in_f) -> (out_f, in_f)
    G = out_f // group_size
    scales_full = np.repeat(scales, group_size, axis=0)   # (out_f, in_f)
    W_dequant = Q * scales_full                          # fp32 matmul, no fp16 materialization for Q
    # NB: Q is materialized as int8 then cast. The real production win is avoiding fp16 W;
    # the int8 cast is unavoidable in pure NumPy. CUDA/NEON versions skip the cast entirely.
    y = x.astype(np.float32) @ W_dequant.T
    return y

This is the honest baseline. The real win in production isn’t avoiding the fp32 cast — it’s that we never hold a fp16 copy of W in memory, which for a 7B model is ~14 GB saved. To get further wins in pure NumPy, chunk along the output dimension:

# quantize/kernels.py (continued)
def fused_dequant_matmul_chunked(Q_packed, scales, shape, x, group_size=128, chunk=256):
    out_f, in_f = shape
    batch = x.shape[0]
    y = np.empty((batch, out_f), dtype=np.float32)
    x32 = x.astype(np.float32, copy=False)
    for i in range(0, out_f, chunk):
        rows = min(chunk, out_f - i)
        Q_chunk = unpack_int4(Q_packed[i:i+rows], in_f).astype(np.float32)
        s_chunk = scales[i // group_size : (i + rows) // group_size]
        s_full = np.repeat(s_chunk, group_size, axis=0)
        W_chunk = Q_chunk * s_full
        y[:, i:i+rows] = x32 @ W_chunk.T
    return y

Chunking keeps the temporary W_chunk small — peak memory is O(chunk · in_f) instead of O(out_f · in_f). This is exactly the pattern used in ExLlamaV2 and the TGI int4 path.

Step 6 — User-facing inference API

# quantize/inference.py
import numpy as np
from .kernels import fused_dequant_matmul_chunked
from .pack import load_quantized

class QuantLinear:
    def __init__(self, path: str, group_size: int = 128, chunk: int = 256):
        self.Q_packed, self.scales, self.shape = load_quantized(path)
        self.group_size = group_size
        self.chunk = chunk

    def forward(self, x: np.ndarray) -> np.ndarray:
        return fused_dequant_matmul_chunked(
            self.Q_packed, self.scales, self.shape, x,
            group_size=self.group_size, chunk=self.chunk,
        )

def load_model(dir_path: str) -> list[QuantLinear]:
    import os
    layers = []
    for fname in sorted(os.listdir(dir_path)):
        if fname.endswith('.npz'):
            layers.append(QuantLinear(os.path.join(dir_path, fname)))
    return layers

Clean, swap-in interface. If you later replace the kernel with a Triton/CUDA version, only kernels.py changes.

Running and Testing It

A real project needs a runnable demo and tests that prove correctness, not just a code dump.

The demo

# demo.py
import numpy as np
from quantize.calibrate import compute_hessian, gptq_quantize
from quantize.pack import pack_int4, save_quantized
from quantize.inference import QuantLinear
import os, tempfile

def main():
    rng = np.random.default_rng(0)
    out_f, in_f = 1024, 1024
    W = rng.normal(0, 0.05, size=(out_f, in_f)).astype(np.float32)
    X = rng.normal(0, 1, size=(256, in_f)).astype(np.float32)

    H = compute_hessian(X)
    Q, scales = gptq_quantize(W, H, group_size=128)
    Q_packed = pack_int4(Q)
    print(f"weight footprint: {Q.nbytes/1e6:.2f} MB (int4)")
    print(f"packed size:      {Q_packed.nbytes/1e6:.2f} MB")

    out_dir = tempfile.mkdtemp()
    path = os.path.join(out_dir, 'layer_0.npz')
    save_quantized(path, Q_packed, scales, Q.shape)

    layer = QuantLinear(path)
    x = rng.normal(0, 1, size=(8, in_f)).astype(np.float16)
    y_int4 = layer.forward(x)
    y_fp32 = x.astype(np.float32) @ W.T
    err = np.abs(y_int4 - y_fp32).mean() / np.abs(y_fp32).mean()
    print(f"relative error:   {err:.4f}")

if __name__ == "__main__":
    main()

Run it:

python demo.py

Expected output: weight footprint ~2.0 MB, packed size ~0.5 MB, relative error somewhere in the 0.005–0.02 range depending on your calibration batch. The 4× memory reduction between fp16 weight matrix and packed int4 is exactly what you want to see.

The tests

# tests/test_quantize.py
import numpy as np
from quantize.calibrate import quantize_group_affine, gptq_quantize, compute_hessian
from quantize.pack import pack_int4, unpack_int4

def test_pack_roundtrip():
    Q = np.array([[-8, -1, 0, 7, -8, 7, 0, 1]], dtype=np.int8)
    packed = pack_int4(Q[:, :8])  # even cols only for simplicity
    out = unpack_int4(packed, 8)
    np.testing.assert_array_equal(out, Q)

def test_group_affine_close_to_fp32():
    rng = np.random.default_rng(42)
    W = rng.normal(0, 0.05, (512, 512)).astype(np.float32)
    Q, scales = quantize_group_affine(W, group_size=128)
    W_hat = Q.astype(np.float32) * scales.repeat(128).reshape(512, -1)
    # SNR should be well above 20 dB for gaussian weights
    signal = (W ** 2).mean()
    noise = ((W - W_hat) ** 2).mean()
    snr_db = 10 * np.log10(signal / noise)
    assert snr_db > 20, f"SNR too low: {snr_db}"

def test_gptq_beats_round_to_nearest():
    rng = np.random.default_rng(7)
    W = rng.normal(0, 0.05, (256, 256)).astype(np.float32)
    # inject outliers
    W[0, 0] = 2.0; W[5, 5] = -3.0
    X = rng.normal(0, 1, (128, 256)).astype(np.float32)
    H = compute_hessian(X)
    Q_rtn, s_rtn = quantize_group_affine(W, 128)
    W_rtn = Q_rtn.astype(np.float32) * s_rtn.repeat(128).reshape(256, -1)
    Q_gptq, s_gptq = gptq_quantize(W, H, 128)
    # approximate reconstruction with group scales
    s_full = s_gptq.repeat(128, axis=0)
    W_gptq = Q_gptq.astype(np.float32) * s_full
    err_rtn = ((W - W_rtn) ** 2).mean()
    err_gptq = ((W - W_gptq) ** 2).mean()
    assert err_gptq <= err_rtn * 1.05, "GPTQ should match or beat RTN"

python -m pytest tests/ -v

If the GPTQ test passes, you’ve demonstrated an empirical fact that shows up in production: second-order calibration matches or beats naive round-to-nearest on the same group size. That’s the result the original GPTQ paper reports.

Extending It: Your Roadmap to Senior-Level

The build above is honest and runnable, but it’s a starting point. Each of the upgrades below turns it into something closer to what a real ML infra team ships. Each is named, scoped, and gives you a concrete PR to point at on your CV.

  1. Triton or CUDA kernel for the fused matmul. Replace the NumPy inner loop with a Triton kernel that does the unpack + scale + matmul in a single pass over global memory, then dispatch by dtype. Why it matters: This is the actual production pattern in vLLM and TGI. Owning one fused kernel is the single strongest signal of low-level ML systems skill.

  2. GPTQ for grouped, asymmetric quantization (int4 with per-group zero-points). Extend pack.py and kernels.py to store and apply per-group zero-points, matching what bitsandbytes and AutoGPTQ do for NF4/int4. Why it matters: Asymmetric quantization squeezes another 10–15% accuracy out of int4 and is what most production engines default to. It also forces you to think carefully about zero-point math in your kernel.

  3. End-to-end pipeline on a real small model. Run your quantizer on facebook/opt-125m or TinyLlama-1.1B, integrate with 🤗 Transformers, and measure perplexity on WikiText-2 against the fp16 baseline. Why it matters: Hiring managers recognize “quantized OPT-125m, perplexity delta 0.05” instantly. It also surfaces real bugs — calibration batch size, attention vs MLP layer handling, embedding layer exclusion — that the toy demo hides.

  4. Benchmarking harness with roofline analysis. Add a bench/ module that times your kernel against torch.matmul on the fp16 reference, reports GB/s achieved, and compares against the theoretical memory-bandwidth roofline (use torch.cuda.get_device_properties for peak bandwidth on GPU, or psutil / /proc/meminfo on CPU). Why it matters: “Kernel hits 78% of HBM peak bandwidth” is a sentence only someone who has actually done the work can write. It also gives you a numbers-driven story for the CV.

  5. Observability hooks: quantization error per layer, SNR histograms, drift detection. Emit per-layer quantization error metrics in a structured format and plot them. Add a CLI flag to dump a side-by-side SNR bar chart. Why it matters: Production quantization is rarely “one shot” — models get re-quantized as data shifts. Showing you thought about monitoring the quantization quality itself is a senior-level signal.

  6. Distributed calibration for models that don’t fit on one device. Use Ray or Dask to shard the Hessian computation across calibration batches on multiple machines, then gather and run GPTQ. Why it matters: This is what quantizing 70B+ models actually requires. Talking credibly about distributed second-order optimization puts you in the conversation for model-serving roles at frontier labs.

Pick two of these for a strong portfolio piece. Pick four and you’re at senior-staff territory for ML infra roles.

Key Takeaways

  • GPTQ is a second-order, column-greedy quantization scheme; the Hessian is what lets it beat naive round-to-nearest on the same group size.
  • Per-group affine scales (typically group_size = 128) are the right granularity trade-off for int4: fine enough to track weight distributions, coarse enough to keep overhead small.
  • The fused dequantize-matmul pattern avoids materializing a fp16 weight matrix at inference time — that 2× memory win over naive dequant-then-matmul is the production argument for kernel fusion.
  • Packing two int4 values per byte with an unsigned offset is the storage format used across llama.cpp, ExLlamaV2, and TGI; it simplifies comparisons and vectorized unpacking.
  • The biggest leap from “toy” to “production-flavored” is replacing the NumPy inner loop with a Triton/CUDA kernel and measuring against a roofline — everything else is engineering polish around that core decision.

Further Reading