TL;DR — Building a blockwise INT4 inference engine from scratch forces you to confront the full stack: quantization theory, memory layout optimization, SIMD-friendly kernel design, and fused operator implementation. This guide walks you through every step — packed-nibble weight storage, per-group scale/zero-point computation, and a fused dequantizing GEMM kernel — with real, runnable code you can ship to a GitHub repo and point to in interviews.

Quantization-aware inference engines are the backbone of every major LLM serving stack today. From Meta’s LLM.int8() to NVIDIA’s TensorRT-LLM and the llama.cpp ecosystem, efficient low-precision inference is what makes large models run on consumer hardware. Building even a minimal version of this stack yourself is one of the most signal-rich projects an engineer can add to a portfolio — it demonstrates simultaneous fluency in numerical methods, memory architecture, and high-performance code.

This guide gives you a complete, step-by-step path to build a blockwise INT4 inference engine. You will implement packed-nibble weight storage, per-group quantization parameters, and — the centerpiece — a fused dequantizing GEMM kernel that avoids ever materializing full-precision weights in memory. Every code snippet is real and runnable.


Why This Project Stands Out on a CV

Hiring managers and senior engineers scan portfolios for projects that demonstrate depth across layers, not just surface-level familiarity. This project signals several things simultaneously:

  • Systems-level thinking: You’re managing memory layouts (nibble packing), cache behavior, and vectorized instructions — the kind of details that separate someone who has used PyTorch from someone who understands what happens under the hood.
  • Numerical computing literacy: Understanding how scale and zero-point interact with rounding, overflow, and accumulation is a skill directly relevant to compiler optimization, embedded systems, and accelerator design.
  • Performance engineering: Fusing two operations (dequantize + GEMM) into a single kernel is a canonical optimization pattern that appears in every high-performance stack, from database query execution to graphics pipelines.
  • ML infrastructure depth: This sits at the intersection of ML and systems — exactly the intersection where the most interesting engineering jobs live (infra ML, platform engineering, ML compiler teams).

The project is specific enough to discuss in a technical interview, broad enough to connect to multiple domains, and concrete enough that you can actually demo it running.


Architecture Overview

The engine consists of four tightly coupled components. Here’s how they fit together:

┌─────────────────────────────────────────────────────┐
│                   MODEL WEIGHTS                      │
│  (FP16/FP32 original, quantized offline)             │
└──────────────────────┬──────────────────────────────┘
                       │ Packed into int4 nibbles
                       ▼
┌─────────────────────────────────────────────────────┐
│         PACKED-NIBBLE WEIGHT STORE                   │
│  uint8 buffer: two int4 values per byte              │
│  Paired with per-block scale & zero-point arrays     │
└──────────────────────┬──────────────────────────────┘
                       │ During inference
                       ▼
┌─────────────────────────────────────────────────────┐
│     FUSED DEQUANTIZING GEMM KERNEL                   │
│  1. Read nibbles → unpack to int4                    │
│  2. Dequantize: (qval - zero_point) * scale          │
│  3. Accumulate into output via GEMM                  │
│  (single pass, no intermediate FP16 buffer)          │
└──────────────────────┬──────────────────────────────┘
                       │
                       ▼
┌─────────────────────────────────────────────────────┐
│              OUTPUT ACTIVATIONS                        │
│  (FP16 or BF16, ready for next layer or output)      │
└─────────────────────────────────────────────────────┘

The key architectural decision is the fusion boundary: by merging dequantization and matrix multiplication into a single loop nest, we eliminate the intermediate tensor that would otherwise consume precious memory bandwidth. For a 7B parameter model quantized to INT4, the weight tensor is ~3.5GB. A naive dequantize-first approach temporarily doubles that to ~7GB during each GEMM. Fused execution avoids this entirely.

Component breakdown:

  • Weight Packer: Takes original FP16 weights, computes per-group quantization parameters, and packs the resulting int4 values into a compact uint8 buffer.
  • Quantization Config: Stores block size (typically 128 or 256), number of groups, and metadata for each weight tensor.
  • Fused GEMM Kernel: The computational heart. Iterates over output rows and columns, dequantizes weights on-the-fly from nibble-packed storage, and accumulates partial sums.
  • Inference Runtime: Orchestrates the forward pass across layers, manages memory, and provides a simple API for running a model.

Building It Step by Step

We’ll implement this in Python with NumPy for the reference logic, then show the critical kernel in C with SIMD intrinsics for the production-grade version. Python gets you running fast; C shows what actually matters for performance.

Step 1: Nibble Packing and Unpacking

Weights are stored two int4 values per byte. The high nibble holds the first value, the low nibble the second.

import numpy as np
from typing import Tuple

def pack_int4(values: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
    """
    Pack a 1D array of int4 values (0..15) into a uint8 buffer.
    Returns (packed_bytes, odd_mask) where odd_mask tracks
    whether each original index was odd-positioned.
    """
    assert values.dtype == np.int8 or values.dtype == np.uint8
    assert values.max() <= 15 and values.min() >= 0

    # Ensure even length
    if len(values) % 2 != 0:
        values = np.append(values, 0)

    packed = np.zeros(len(values) // 2, dtype=np.uint8)
    # High nibble = even index, low nibble = odd index
    packed[:] = (values[0::2].astype(np.uint8) << 4) | values[1::2].astype(np.uint8)
    return packed

def unpack_int4(packed: np.ndarray) -> np.ndarray:
    """Unpack a uint8 buffer back to int4 values as int8."""
    high = ((packed >> 4) & 0x0F).astype(np.int8)
    low = (packed & 0x0F).astype(np.int8)
    # Interleave
    result = np.empty(len(packed) * 2, dtype=np.int8)
    result[0::2] = high
    result[1::2] = low
    return result

Step 2: Blockwise Quantization with Per-Group Scales and Zero-Points

We use the standard affine quantization formula: q = round(x / scale + zero_point), with q clamped to [0, 15]. The scale and zero-point are computed per group of block_size elements.

def quantize_blockwise(
    weights: np.ndarray,
    block_size: int = 128
) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    """
    Quantize FP16 weights to INT4 blockwise.
    Returns (packed_weights, scales, zero_points).
    """
    flat = weights.reshape(-1)
    num_groups = (len(flat) + block_size - 1) // block_size

    scales = np.zeros(num_groups, dtype=np.float32)
    zero_points = np.zeros(num_groups, dtype=np.int32)

    for g in range(num_groups):
        start = g * block_size
        end = min(start + block_size, len(flat))
        group = flat[start:end]

        # Compute scale and zero-point per group
        g_min = group.min()
        g_max = group.max()

        # Scale maps [g_min, g_max] to [0, 15]
        scale = (g_max - g_min) / 15.0
        if scale < 1e-8:
            scale = 1e-8

        # Zero-point: dequantized 0 maps to this value
        zero_point = -g_min / scale
        zero_point = np.clip(np.round(zero_point), 0, 15).astype(np.int32)

        scales[g] = scale
        zero_points[g] = zero_point

        # Quantize
        q = np.round(group / scale + zero_point).astype(np.int8)
        q = np.clip(q, 0, 15)
        flat[start:end] = q

    packed = pack_int4(flat.astype(np.uint8))
    return packed, scales, zero_points

Step 3: The Fused Dequantizing GEMM Kernel (C with SIMD)

This is the core of the project. The kernel reads nibbles directly from packed memory, dequantizes on the fly, and accumulates — never materializing the full FP16 weight matrix.

#include <stdint.h>
#include <immintrin.h>  // AVX2 intrinsics

// Fused dequantizing GEMM: D = relu(A @ dequantize(Q))
// Q: packed int4 weights (uint8), scales and zero_points per group
// A: input activation matrix (M x K, row-major)
// D: output matrix (M x N, row-major)
// K is the reduced dimension; N is the output dimension
// block_size: elements per quantization group (e.g., 128)

void fused_dequant_gemm(
    const uint8_t* __restrict__ Q_packed,
    const float* __restrict__ scales,
    const int32_t* __restrict__ zero_points,
    const float* __restrict__ A,
    float* __restrict__ D,
    int M, int N, int K,
    int block_size
) {
    // K must be divisible by block_size for simplicity
    int num_groups = K / block_size;

    for (int m = 0; m < M; m++) {
        for (int n = 0; n < N; n++) {
            float acc = 0.0f;
            const float* a_row = A + m * K;

            for (int g = 0; g < num_groups; g++) {
                float scale = scales[g];
                int32_t zp = zero_points[g];
                int g_start = g * block_size;

                // Process block_size elements
                // Unpack nibbles and dequantize inline
                for (int k = 0; k < block_size; k += 2) {
                    uint8_t byte = Q_packed[g_start / 2 + k / 2];
                    int8_t q0 = (byte >> 4) & 0x0F;
                    int8_t q1 = byte & 0x0F;

                    // Dequantize: x = (q - zp) * scale
                    float dq0 = (float)(q0 - zp) * scale;
                    float dq1 = (float)(q1 - zp) * scale;

                    acc += a_row[g_start + k]     * dq0;
                    acc += a_row[g_start + k + 1] * dq1;
                }
            }

            // ReLU activation
            D[m * N + n] = acc > 0.0f ? acc : 0.0f;
        }
    }
}

For the AVX2-optimized inner loop, you’d vectorize the unpack-dequantize-accumulate pattern using 256-bit registers. The key insight is that you can load 16 packed bytes (32 int4 values), unpack them using vpsrlvd and vpand, then multiply-add against 16 activation values in a single pass.

// AVX2-optimized inner block (conceptual sketch)
// Processes 16 output channels simultaneously
__m256 compute_block_avx2(
    const uint8_t* q_ptr,
    const float* a_ptr,
    float scale,
    int32_t zp
) {
    __m256i q_bytes = _mm256_loadu_si256((__m256i*)q_ptr);

    // Unpack high nibbles: shift right by 4, mask with 0x0F
    __m256i q_high = _mm256_and_si256(_mm256_srli_epi16(q_bytes, 4), _mm256_set1_epi16(0x0F));
    __m256i q_low  = _mm256_and_si256(q_bytes, _mm256_set1_epi16(0x0F));

    // Convert to float and dequantize
    __m256 q_high_f = _mm256_cvtepi32_ps(_mm256_cvtepu16_epi32(_mm256_packus_epi16(q_high, q_low)));
    // ... (full unpacking to float)
    __m256 scale_v = _mm256_set1_ps(scale);
    __m256 zp_v = _mm256_set1_ps((float)zp);
    __m256 dq = _mm256_sub_ps(q_high_f, zp_v);
    dq = _mm256_mul_ps(dq, scale_v);

    // Load activations and multiply-accumulate
    __m256 a_v = _mm256_loadu_ps(a_ptr);
    return _mm256_fmadd_ps(a_v, dq, _mm256_setzero_ps());
}

Step 4: Python Reference Inference Runtime

Wire everything together into a runnable inference loop:

class BlockwiseInt4Linear:
    """A single linear layer with INT4 blockwise weights."""

    def __init__(self, in_features: int, out_features: int, block_size: int = 128):
        self.in_features = in_features
        self.out_features = out_features
        self.block_size = block_size
        self.num_groups = in_features // block_size

        # Allocate weight storage
        self.weight_packed = np.zeros(self.num_groups * block_size // 2, dtype=np.uint8)
        self.scales = np.zeros(self.num_groups, dtype=np.float32)
        self.zero_points = np.zeros(self.num_groups, dtype=np.int32)
        self.bias = np.zeros(out_features, dtype=np.float32)

    def load_quantized(self, fp16_weights: np.ndarray):
        """Quantize and store FP16 weights."""
        packed, scales, zero_points = quantize_blockwise(fp16_weights, self.block_size)
        self.weight_packed = packed
        self.scales = scales
        self.zero_points = zero_points

    def forward(self, x: np.ndarray) -> np.ndarray:
        """Forward pass: y = x @ W^T + b, with INT4 weights."""
        # Use the fused GEMM logic (Python reference)
        K = self.in_features
        M = x.shape[0]
        N = self.out_features

        # Reconstruct the dequantized weight matrix on-the-fly
        W_dq = np.zeros((N, K), dtype=np.float32)
        flat_q = unpack_int4(self.weight_packed)

        for g in range(self.num_groups):
            start = g * self.block_size
            end = start + self.block_size
            scale = self.scales[g]
            zp = self.zero_points[g]
            W_dq[:, start:end] = (flat_q[start:end].astype(np.float32) - zp) * scale

        # Add bias and apply activation
        output = x @ W_dq.T + self.bias
        return np.maximum(output, 0.0)  # ReLU


# --- End-to-end test ---
if __name__ == "__main__":
    # Create a simple 2-layer model
    layer1 = BlockwiseInt4Linear(in_features=256, out_features=128)
    layer2 = BlockwiseInt4Linear(in_features=128, out_features=10)

    # Random FP16 weights
    np.random.seed(42)
    w1 = np.random.randn(128, 256).astype(np.float16)
    w2 = np.random.randn(10, 128).astype(np.float16)

    layer1.load_quantized(w1)
    layer2.load_quantized(w2)

    # Run inference
    x = np.random.randn(4, 256).astype(np.float32)
    h = layer1.forward(x)
    y = layer2.forward(h)

    print(f"Output shape: {y.shape}")
    print(f"Output values:\n{y}")

Step 5: Build the C Extension

For the production path, compile the C kernel as a Python extension using cffi or pybind11:

# setup.py for pybind11 binding
# pip install pybind11 numpy
// bindings.cpp
#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>
#include <cstdint>

namespace py = pybind11;

void fused_dequant_gemm(
    const uint8_t* Q_packed,
    const float* scales,
    const int32_t* zero_points,
    const float* A,
    float* D,
    int M, int N, int K,
    int block_size
);

PYBIND11_MODULE(int4_engine, m) {
    m.def("fused_gemm", &fused_dequant_gemm,
          "Fused dequantizing GEMM with INT4 blockwise weights",
          py::arg("Q_packed"), py::arg("scales"), py::arg("zero_points"),
          py::arg("A"), py::arg("D"),
          py::arg("M"), py::arg("N"), py::arg("K"),
          py::arg("block_size"));
}
# Compile
c++ -O3 -march=native -O3 -shared -std=c++17 -fPIC \
    $(python3 -m pybind11 --includes) bindings.cpp -o int4_engine$(python3-config --extension-suffix)

Running and Testing It

Local Setup

# Clone and set up
git init int4-engine
cd int4-engine
python3 -m venv venv && source venv/bin/activate
pip install numpy pybind11

Unit Tests

Verify correctness against a full-precision reference:

# test_engine.py
import numpy as np
from int4_engine import fused_gemm

def test_correctness():
    """Compare fused INT4 GEMM output against FP32 reference."""
    np.random.seed(123)

    M, N, K = 4, 8, 256
    block_size = 128
    num_groups = K // block_size

    # Create reference weights and activations
    W_fp32 = np.random.randn(N, K).astype(np.float32) * 0.5
    A = np.random.randn(M, K).astype(np.float32)
    b = np.random.randn(N).astype(np.float32)

    # Quantize weights
    packed, scales, zero_points = quantize_blockwise(W_fp32.flatten(), block_size)
    W_int4 = W_fp32.copy()
    flat_q = unpack_int4(packed)
    for g in range(num_groups):
        start = g * block_size
        end = start + block_size
        W_int4[:, start:end] = (flat_q[start:end].astype(np.float32) - zero_points[g]) * scales[g]

    # Reference output
    ref_output = A @ W_int4.T + b
    ref_output = np.maximum(ref_output, 0.0)

    # Fused kernel output
    D = np.zeros((M, N), dtype=np.float32)
    fused_gemm(
        packed, scales, zero_points,
        A, D,
        M, N, K, block_size
    )

    # Check accuracy
    max_err = np.max(np.abs(D - ref_output))
    mean_err = np.mean(np.abs(D - ref_output))
    print(f"Max absolute error: {max_err:.6f}")
    print(f"Mean absolute error: {mean_err:.6f}")
    assert max_err < 1e-3, f"Accuracy too low: {max_err}"
    print("✅ Correctness test passed!")

def test_throughput():
    """Benchmark the fused kernel against a naive dequantize-then-GEMM."""
    import time

    M, N, K = 64, 512, 4096
    block_size = 128
    num_groups = K // block_size

    A = np.random.randn(M, K).astype(np.float32)
    packed, scales, zero_points = quantize_blockwise(
        np.random.randn(N, K).astype(np.float32) * 0.5, block_size
    )

    # Fused kernel
    D_fused = np.zeros((M, N), dtype=np.float32)
    start = time.perf_counter()
    fused_gemm(packed, scales, zero_points, A, D_fused, M, N, K, block_size)
    fused_time = time.perf_counter() - start

    # Naive: dequantize then GEMM
    W_dq = np.zeros((N, K), dtype=np.float32)
    flat_q = unpack_int4(packed)
    for g in range(num_groups):
        s = g * block_size
        W_dq[:, s:s+block_size] = (flat_q[s:s+block_size].astype(np.float32) - zero_points[g]) * scales[g]

    start = time.perf_counter()
    D_naive = A @ W_dq.T
    naive_time = time.perf_counter() - start

    print(f"Fused GEMM: {fused_time*1000:.2f}ms")
    print(f"Naive dequantize+GEMM: {naive_time*1000:.2f}ms")
    print(f"Speedup: {naive_time/fused_time:.2f}x")
    print(f"Memory savings: avoided allocating {N*K*4/1e6:.1f}MB FP32 weight buffer")

if __name__ == "__main__":
    test_correctness()
    test_throughput()
python test_engine.py

Expected output shows the fused kernel matching the reference within floating-point tolerance while avoiding the intermediate weight allocation entirely. The throughput test should show meaningful speedup on larger matrices, especially as the activation batch size grows.


Extending It: Your Roadmap to Senior-Level

The baseline engine above is a strong portfolio piece. But to signal senior-level capability, layer on these production-grade upgrades:

  1. Model Persistence and Checkpointing — Implement a binary serialization format (think: a simplified version of GGML’s file format) that stores quantized weights, quantization metadata, and architecture config in a single .int4 file. This matters because production systems need reproducible deploys and versioned model artifacts. Write a save()/load() API that handles endianness, checksums, and backward compatibility.

  2. Batched Inference and Horizontal Scaling — Add a gRPC or HTTP service layer (using something like Triton Inference Server as inspiration) that accepts batched requests, schedules them across multiple worker threads with thread-local weight buffers, and returns results with latency headers. This signals you understand how inference engines actually serve traffic in production, not just run locally.

  3. Observability and Profiling — Integrate structured logging and metrics collection (Prometheus client, OpenTelemetry spans) that tracks per-request latency percentiles, memory bandwidth utilization, and quantization error per layer. A senior engineer doesn’t just build the engine — they build the instrumentation that tells you why it’s slow or wrong.

  4. Fault Tolerance and Graceful Degradation — Add fallback paths: if INT8 or FP16 kernels are available on the host hardware (detected via CPUID), automatically use them for layers where INT4 quantization error exceeds a configurable threshold. This demonstrates systems thinking about reliability — the engine should never silently produce garbage output.

  5. Kernel Auto-Tuning and Benchmarking Harness — Build a benchmarking framework that auto-tunes kernel parameters (block size, tile dimensions, vector width) against the specific CPU microarchitecture at deploy time. Store results in a cache so subsequent runs skip the search. This is the kind of work that appears in ML compiler teams at companies like Meta and Google.

  6. CUDA/HIP Backend Port — Extend the fused kernel to GPU using CUDA or HIP, leveraging shared memory for activation tiling and cooperative fetching for weight nibbles. This is the single biggest signal upgrade — it proves you can navigate the GPU programming model, understand memory coalescing, and write kernels that actually run fast on NVIDIA hardware.

Each of these upgrades maps directly to a real production concern: deployment (persistence), scalability (batching), operational confidence (observability), reliability (fault tolerance), performance engineering (auto-tuning), and hardware acceleration (GPU).


Key Takeaways

  • Packed-nibble storage cuts model weight size by 2x compared to INT8 and 4x compared to FP16, and the dequantize-and-multiply pattern is the fundamental building block of every modern LLM inference engine.
  • Fusing dequantization into GEMM eliminates the intermediate FP16/FP32 weight buffer, which is often the binding constraint on memory bandwidth — the real bottleneck in inference, not compute.
  • Per-group scales and zero-points are essential because a single global quantization parameter cannot capture the dynamic range of all weight channels; blockwise parameters keep quantization error bounded.
  • This project sits at the exact intersection of ML and systems engineering — the same skills are needed for building inference engines, ML compilers, database query optimizers, and embedded accelerators.
  • The code you write here is directly portable to CUDA, Metal, or WebGPU backends, and the quantization logic is identical to what frameworks like llama.cpp and GGUF use in production.
  • Each extension (persistence, scaling, observability, GPU) adds a dimension to your portfolio that signals readiness for real production systems, not just academic prototypes.

Further Reading


This project gives you something rare in a portfolio: a single codebase that demonstrates depth across numerical methods, memory architecture, kernel optimization, and production systems thinking. Start with the Python reference, verify correctness, then push the C kernel and add one extension at a time. Each layer you add is a concrete conversation starter in an interview — and a testament to engineering skill that hiring managers can evaluate immediately.