Lesson 26: Quantization and Fixed-Point Inference

Every parameter of the model you trained in 23-putting-it-all-together is a 64-bit float. A 70-billion-parameter model stored the same way at fp16 occupies 140 GB, and during generation every one of those bytes is read from memory once per emitted token — so the number of bits per weight is a throughput and memory knob as real as \(d_{\text{model}}\). This lesson re-trains the capstone model, quantises its weights to \(b \in \{2, \dots, 8\}\) bits with rustlab's integer and fixed-point types, and measures what each bit costs: in signal-to-noise ratio per tensor, in perplexity, and in the KL divergence of the next-token distribution. Read together, those measurements are the rate–distortion curve of a language model — the last row of the toolkit table in 00-the-llm-as-a-system.

Learning Objectives

• Implement symmetric uniform quantisation of a weight tensor — scale, round, clip, dequantise — with per-tensor and per-row scales, and read it as a Q-format whose step is tuned to the data.
• Predict the SNR of a quantised tensor from the additive uniform-noise model (\(\Delta^2/12\), 6.02 dB per bit, minus a crest-factor penalty) and measure where the model holds and where it breaks.
• Measure the rate–distortion curve of the capstone model: perplexity and \(\mathrm{KL}(p \,\|\, p_q)\) versus bits per weight, and locate its knee.
• Explain why activations are harder to quantise than weights, why saturation beats wraparound, and how a quantisation-aware training loop differs from post-training quantisation.
• Compute the KV-cache memory and the bandwidth-bound decode rate of a 70B model at 16, 8 and 4 bits.

Background

The trained model, its 32-token BPE corpus and its 300 parameters from 23-putting-it-all-together; perplexity and cross-entropy in nats and bits from 20-perplexity-and-evaluation; KL divergence from 03-cross-entropy-loss; the KV cache from 21-sampling-and-generation and GQA from 24-modern-architectural-variants. From engineering: the uniform quantiser of an ADC, two's-complement integers, and Q-format fixed-point notation. Loss is computed in nats and reported in bits (\(\div \ln 2\)) where labelled.

The Model Under Test

Theory

The capstone's single-block transformer has 300 parameters: 272 of them sit in eight weight matrices — \(\mathbf{E}\) (\(18 \times 4\)), \(\mathbf{W}_Q, \mathbf{W}_K, \mathbf{W}_V, \mathbf{W}_O\) (\(4 \times 4\)), \(\mathbf{W}_1\) (\(4 \times 8\)), \(\mathbf{W}_2\) (\(8 \times 4\)), \(\mathbf{W}_U\) (\(4 \times 18\)) — and the remaining 28 are LayerNorm affines and FFN biases. Those eight matrices are what gets quantised in this lesson; the 28 vectors stay in full precision, as they do in every production quantiser (in a real model they are well under 0.1 % of the parameters). The hidden setup block above reproduces the capstone's training run exactly — same seed, same 600 AdamW steps — so the weights below are the weights Lesson 23 ended with. It also defines mat_k(P, k), a hidden helper that returns the \(k\)-th of the eight matrices in the order above (rustlab has no struct-field indexing), which the visible blocks use to loop over them.

Example — Reproduce the capstone

[lo_fp, L_fp, cache_fp] = transformer_forward(tokens, mask_train, P);
print("re-trained capstone: params =", n_params, "  L =", L_fp, "  PPL =", exp(L_fp));
n_q = 0;
ranges = zeros(8);
for k = 1:8
  W = mat_k(P, k);
  n_q = n_q + numel(W);
  ranges(k) = max(max(abs(W)));
end
print("quantised parameters:", n_q, "  max|W| per matrix:", ranges);

re-trained capstone: params = 300   L = 0.00008169636170878258   PPL = 1.0000816996989474
quantised parameters: 272   max|W| per matrix: [1×8]  2.008726  1.961393  1.477834  0.722783  0.943531  2.144984  1.989670  1.578401

\(\mathcal{L} = 8.17 \times 10^{-5}\) nats and \(\mathrm{PPL} = 1.00008\) — the capstone's numbers. The weights started at \(\mathcal{N}(0, 0.3^2)\) and training stretched them to ranges between \(0.72\) and \(2.14\); that spread is the first thing a quantiser has to cope with.

Why Quantise

Theory

Two costs scale with bits per weight. Memory: a model of \(N\) parameters at \(b\) bits occupies \(N b / 8\) bytes — 140 GB for \(N = 7 \times 10^{10}\) at fp16, more than most single accelerators hold. Bandwidth: autoregressive decode does one matrix–vector product per layer per token, so every weight is read from memory exactly once per generated token and almost no arithmetic is done on it. The token rate is therefore bounded by

\[\text{tok/s} \;\le\; \frac{\text{memory bandwidth}}{\text{bytes read per token}} = \frac{B}{N b / 8 + \text{KV-cache bytes}},\]

an Exact upper bound (arithmetic is assumed free). 24-modern-architectural-variants applies the same bound to the KV cache to motivate GQA; here the weight term dominates. Halving \(b\) doubles the bound — which is why bits per weight is a design parameter, not an implementation detail.

Example — The byte budget of a 70B model

N_llm = 70e9;                    % parameters
BW = 3.35e12;                    % bytes/s — one H100 SXM (HBM3)
w_bits = [16, 8, 4];
w_GB = N_llm * w_bits / 8 / 1e9;
tok_bound = BW ./ (w_GB * 1e9);
print("bits per weight:", w_bits);
print("weights (GB):   ", w_GB);
print("tok/s bound:    ", tok_bound);

bits per weight: [1×3]  16.000000  8.000000  4.000000
weights (GB):    [1×3]  140.000000  70.000000  35.000000
tok/s bound:     [1×3]  23.928571  47.857143  95.714286

At fp16 the weights alone cap decode at \(23.9\) tokens per second on one device; int4 lifts the same bound to \(95.7\). Whether the quality survives the trip from 16 bits to 4 is the question the rest of the lesson measures.

Uniform Quantisation of a Tensor

Theory

A symmetric uniform quantiser with \(b\) bits maps a real tensor \(\mathbf{W}\) to integers \(\mathbf{q}\) in \([-q_{\max}, q_{\max}]\), \(q_{\max} = 2^{b-1} - 1\), through one real scale \(s\):

\[s = \frac{\max \lvert \mathbf{W} \rvert}{q_{\max}}, \qquad \mathbf{q} = \operatorname{clip}\!\left(\operatorname{round}(\mathbf{W} / s),\, -q_{\max},\, q_{\max}\right), \qquad \hat{\mathbf{W}} = s\,\mathbf{q}.\]

The step \(\Delta = s\) is the resolution; \(\hat{\mathbf{W}}\) is the dequantised tensor that the forward pass actually multiplies with. Memory holds \(\mathbf{q}\) at \(b\) bits per entry plus the scale; the matmul reads \(\mathbf{q}\), multiplies by \(s\), and proceeds as before — nothing else in the model changes.

Choosing \(s\) from \(\max\lvert\mathbf{W}\rvert\) guarantees nothing saturates, at the price that one outlier sets the step for the whole tensor. The remedy is finer scale granularity: a scale per row of \(\mathbf{W}\) (for \(\mathbf{E}\) a row is one token's embedding; for a projection it is one input channel) so that a small-range row keeps a small step. The cost is one extra float per row.

Fixed-point arithmetic is the special case \(s = 2^{-f}\): a Q-format with \(f\) fractional bits. rustlab's qfmt(word, frac) and quantize implement exactly that grid, int8(x) rounds half away from zero and saturates by default, and snr(x, x_q) returns \(10 \log_{10}\) of signal power over error power in dB.

Example — int8 and the Q-format view of the LM head

W = P.W_U;
n_WU = numel(W);
s = max(max(abs(W))) / 127;
q = int8(W / s);                                    % rounds half away from zero, saturates at +-127
W_hat = double(q) * s;
print("max|W_U| =", max(max(abs(W))), "  step s =", s, "  class(q) =", class(q));
print("W_U(1, 1:6) =", W(1, 1:6));
print("q(1, 1:6)   =", q(1, 1:6));
print("entries at +-127:", sum(sum(abs(double(q)) == 127)), "   int8 SNR (dB):", snr(reshape(W, 1, n_WU), reshape(W_hat, 1, n_WU)));
fmt_07 = qfmt(8, 7, "round", "saturate");           % eight bits, seven fractional: range [-1, 1)
fmt_16 = qfmt(8, 6, "round", "saturate");           % one integer bit: range [-2, 2), step 2^-6
snr_07 = snr(reshape(W, 1, n_WU), reshape(quantize(W, fmt_07), 1, n_WU));
snr_16 = snr(reshape(W, 1, n_WU), reshape(quantize(W, fmt_16), 1, n_WU));
print(fmt_07, " saturated entries:", sum(sum(abs(W) >= 1)), "  SNR", snr_07, "dB");
print(fmt_16, " saturated entries:", sum(sum(abs(W) >= 2)), "  SNR", snr_16, "dB   (step", 2 ^ -6, ")");

max|W_U| = 1.5784014486102542   step s = 0.012428357863072868   class(q) = int8
W_U(1, 1:6) = [1×6]  0.433060  0.447733  0.526826  0.709158  0.512755  0.419608
q(1, 1:6)   = [1×6]  35.000000  36.000000  42.000000  57.000000  41.000000  34.000000
entries at +-127: 1    int8 SNR (dB): 45.94462791186676
QFmt<8-bit Q0.7, round=round, overflow=saturate>  saturated entries: 7   SNR 18.639546316414158 dB
QFmt<8-bit Q1.6, round=round, overflow=saturate>  saturated entries: 0   SNR 43.34727862310058 dB   (step 0.015625 )

The step is \(s = 0.0124\); the entry \(W_U(1,1) = 0.4331\) becomes the integer \(35\) and dequantises to \(0.4350\) — an error of \(0.00193\), less than half a step. Exactly one entry (the maximum) sits at the rail. The fixed-point view makes the trade-off explicit. qfmt(8, 7) — rustlab prints it Q0.7; texts that count the sign bit call it Q1.7 — cannot represent \(1.58\), so \(7\) entries saturate and the SNR collapses to \(18.6\) dB. Giving up one fractional bit (Q1.6) fits the range and recovers \(43.3\) dB. The per-tensor scale above does better still (step \(0.0124\) instead of \(2^{-6} = 0.0156\)) because it uses all 127 codes for the actual range rather than the next power of two — a per-tensor scale is a Q-format whose step is tuned to the data, which is the entire difference between classical fixed-point DSP and "integer quantisation" in machine learning.

Example — Per-tensor versus per-row scales

The general quantiser for any \(b\), and its per-row variant. snr_mat flattens a matrix because snr takes vectors.

function [Wq, s] = quant_tensor(W, b)
  q_max = 2 ^ (b - 1) - 1;
  s = max(max(abs(W))) / q_max;
  q = min(max(round(W / s), -q_max), q_max);
  Wq = q * s;
end
function [Wq, s] = quant_rows(W, b)
  q_max = 2 ^ (b - 1) - 1;
  s = max(abs(W), [], 2) / q_max;                % one scale per row (a column vector)
  q = min(max(round(W ./ s), -q_max), q_max);
  Wq = q .* s;
end
function d = snr_mat(W, Wq)
  n = numel(W);
  d = snr(reshape(W, 1, n), reshape(Wq, 1, n));
end

b_demo = 4;
snr_t4 = zeros(8);  snr_r4 = zeros(8);
for k = 1:8
  W = mat_k(P, k);
  [Wt, st] = quant_tensor(W, b_demo);  snr_t4(k) = snr_mat(W, Wt);
  [Wr, sr] = quant_rows(W, b_demo);    snr_r4(k) = snr_mat(W, Wr);
end
row_max_E = max(abs(P.E), [], 2);
print("4-bit SNR per-tensor (dB):", snr_t4);
print("4-bit SNR per-row    (dB):", snr_r4);
figure();
subplot(1, 2, 1)
bar(1:vocab, row_max_E')
title("Row ranges of E: max |E(i, :)| per token")
xlabel("token id i"); ylabel("max |E(i, :)|")
subplot(1, 2, 2)
bar(1:8, [snr_t4', snr_r4'])
legend("per-tensor", "per-row")
ylim([0, 30])
title("4-bit SNR per weight matrix (1 = E, 2-5 = Wq Wk Wv Wo, 6-7 = W1 W2, 8 = W_U)")
xlabel("matrix index"); ylabel("SNR (dB)")

4-bit SNR per-tensor (dB): [1×8]  17.959081  22.175934  22.639553  22.930464  22.891831  21.046100  21.888738  20.832276
4-bit SNR per-row    (dB): [1×8]  25.941052  27.473702  26.253828  26.628783  24.780850  22.369411  24.097818  23.107567

Quantisation Noise

Theory

Write the quantised tensor as signal plus error, \(\hat{\mathbf{W}} = \mathbf{W} + \mathbf{e}\). When the step \(\Delta\) is small compared with the spread of the data, the rounding error is well modelled as independent of the signal and uniform on \([-\Delta/2, \Delta/2]\), with zero mean and variance

\[\sigma_e^2 = \frac{1}{\Delta}\int_{-\Delta/2}^{\Delta/2} e^2 \, de = \frac{\Delta^2}{12}.\]

With \(\Delta = \max\lvert\mathbf{W}\rvert / (2^{b-1} - 1) \approx 2^{1-b} \max\lvert\mathbf{W}\rvert\) and signal power \(P_W = \operatorname{mean}(W_{ij}^2)\),

\[\mathrm{SNR}_{\text{dB}} = 10\log_{10}\frac{P_W}{\Delta^2 / 12} \approx 6.02\,b + 4.77 - 20\log_{10} c, \qquad c = \frac{\max\lvert\mathbf{W}\rvert}{\sqrt{P_W}},\]

where \(c\) is the crest factor. Every extra bit halves \(\Delta\) and buys 6.02 dB — the ADC rule, Exact under the uniform-noise model — and a heavy-tailed tensor (large \(c\)) pays a fixed penalty at every \(b\): the outlier sets the step, the bulk of the entries sit far below full scale. The model breaks at small \(b\), where the error is no longer independent of the signal (at \(b = 2\) the three codes \(\{-1, 0, 1\}\) carry almost none of the data), and it slightly under-predicts the slope because \(2^{b-1} - 1\) shrinks the step a little faster than \(2\times\) per bit at low \(b\).

Example — SNR versus bits for every weight matrix

bits = 2:8;
n_b = length(bits);
SNR_t = zeros(n_b, 8);
SNR_r = zeros(n_b, 8);
for i = 1:n_b
  for k = 1:8
    W = mat_k(P, k);
    [Wt, st] = quant_tensor(W, bits(i));  SNR_t(i, k) = snr_mat(W, Wt);
    [Wr, sr] = quant_rows(W, bits(i));    SNR_r(i, k) = snr_mat(W, Wr);
  end
end
crest_E = max(max(abs(P.E))) / sqrt(mean(mean(P.E .^ 2)));
model_E = 6.02 * bits + 4.77 - 20 * log10(crest_E);
b_fine = linspace(2, 8, 61);                                   % fine grid so the dashed model line draws cleanly
gain_per_bit = SNR_t(2:n_b, :) - SNR_t(1:(n_b - 1), :);        % dB gained by each added bit
mean_gain = zeros(n_b - 1);
for i = 1:(n_b - 1)
  mean_gain(i) = mean(gain_per_bit(i, :));
end
print("crest factor of E:", crest_E, "  model intercept (dB):", 4.77 - 20 * log10(crest_E));
print("E: measured SNR", reshape(SNR_t(:, 1), 1, n_b), "  model", model_E);
print("mean dB per added bit (3..8):", mean_gain);
figure();
subplot(1, 2, 1)
hold("on")
for k = 1:8
  plot(bits, SNR_t(:, k)', "color", "gray")
end
plot(bits, SNR_t(:, 1)', "color", "blue", "label", "E measured")
plot(b_fine, 6.02 * b_fine + 4.77 - 20 * log10(crest_E), "color", "red", "style", "dashed", "label", "E model: 6.02 b + 4.77 - 20 log10(c)")
hold("off")
title("Per-tensor SNR vs bits (gray: the other seven matrices)")
xlabel("bits per weight b"); ylabel("SNR (dB)")
subplot(1, 2, 2)
plot(bits(2:n_b), mean_gain, "color", "purple", "label", "mean over 8 matrices")
hold("on")
yline(6.02, "red", "6.02 dB per bit")
hold("off")
title("dB gained by the b-th bit")
xlabel("bits per weight b"); ylabel("SNR(b) - SNR(b-1)  (dB)")

crest factor of E: 2.8769919967334148   model intercept (dB): -4.408773075280628
E: measured SNR [1×7]  3.901893  11.132618  17.959081  25.080645  31.795283  37.245509  42.976012   model [1×7]  7.631227  13.651227  19.671227  25.691227  31.711227  37.731227  43.751227
mean dB per added bit (3..8): [1×6]  9.898346  7.277465  6.469852  6.664511  5.572157  6.544953

Rate–Distortion for a Language Model

Theory

SNR measures distortion in the weights; what matters is distortion in the output. Run the quantised model — dequantised \(\hat{\mathbf{W}}\) in the unchanged transformer_forward — on the corpus and measure two things against the fp64 model: the loss (or perplexity), and the KL divergence between the two next-token distributions averaged over positions,

\[D(b) = \frac{1}{T}\sum_{t=1}^{T} \mathrm{KL}\!\left(p_t \,\|\, p^{(b)}_t\right) \quad \text{bits}, \qquad p_t = \operatorname{softmax}(\mathbf{z}_t), \; p^{(b)}_t = \operatorname{softmax}(\hat{\mathbf{z}}_t).\]

KL is the right distortion for a language model because it compares the quantised model with the model, not with the data — on a real corpus a quantised model can match the fp perplexity to three digits while assigning visibly different probabilities. (On this memorised corpus \(p_t\) is essentially one-hot, so \(D(b)\) nearly equals the quantised loss in bits; the two diverge on any corpus with real entropy.) The rate is bits per weight; the total is \(272\,b\) bits for the eight matrices plus the scales. Plotting distortion against rate gives the rate–distortion curve, and its knee — the smallest \(b\) beyond which distortion stops falling meaningfully — is the operating point a deployment picks.

Example — Perplexity and KL versus bits per weight

quant_params applies one quantiser to all eight matrices (a hidden helper — one assignment per matrix, because rustlab cannot loop over struct fields) by calling quant_any(W, b, per_row), a second hidden helper that dispatches to quant_rows or quant_tensor; mean_kl averages kl_bits over the \(T - 1\) predicting positions (the last logit row predicts nothing).

L_t = zeros(n_b);  L_r = zeros(n_b);  KL_t = zeros(n_b);  KL_r = zeros(n_b);
for i = 1:n_b
  Pq = quant_params(P, bits(i), 0);
  [lo_q, L_q, c] = transformer_forward(tokens, mask_train, Pq);
  L_t(i) = L_q;  KL_t(i) = mean_kl(lo_fp, lo_q, T - 1);
  Pq = quant_params(P, bits(i), 1);
  [lo_q, L_q, c] = transformer_forward(tokens, mask_train, Pq);
  L_r(i) = L_q;  KL_r(i) = mean_kl(lo_fp, lo_q, T - 1);
  print(sprintf("%d bits  PPL per-tensor %11.5f  per-row %9.5f   KL(bits) per-tensor %9.2e  per-row %9.2e", bits(i), exp(L_t(i)), exp(L_r(i)), KL_t(i), KL_r(i)));
end
PPL_t = exp(L_t);  PPL_r = exp(L_r);
figure();
subplot(1, 2, 1)
semilogy(bits, PPL_t, "color", "blue", "label", "per-tensor")
hold("on")
semilogy(bits, PPL_r, "color", "red", "label", "per-row")
yline(exp(L_fp), "gray", "fp64 model")
hold("off")
title("Perplexity vs bits per weight")
xlabel("bits per weight b"); ylabel("PPL (log scale)")
subplot(1, 2, 2)
semilogy(bits, KL_t, "color", "blue", "label", "per-tensor")
hold("on")
semilogy(bits, KL_r, "color", "red", "label", "per-row")
hold("off")
title("Mean KL(p || p_q) vs bits per weight")
xlabel("bits per weight b"); ylabel("KL (bits, log scale)")

2 bits  PPL per-tensor 11177.06975  per-row 167.77026   KL(bits) per-tensor  1.34e+01  per-row  7.39e+00
3 bits  PPL per-tensor    14.72530  per-row   1.65351   KL(bits) per-tensor  3.88e+00  per-row  7.25e-01
4 bits  PPL per-tensor     3.03257  per-row   1.11390   KL(bits) per-tensor  1.60e+00  per-row  1.55e-01
5 bits  PPL per-tensor     1.06840  per-row   1.00272   KL(bits) per-tensor  9.53e-02  per-row  3.78e-03
6 bits  PPL per-tensor     1.00023  per-row   1.00009   KL(bits) per-tensor  2.13e-04  per-row  1.62e-05
7 bits  PPL per-tensor     1.00012  per-row   1.00011   KL(bits) per-tensor  6.36e-05  per-row  3.62e-05
8 bits  PPL per-tensor     1.00008  per-row   1.00008   KL(bits) per-tensor  4.55e-06  per-row  1.74e-06

The knee is sharp. With per-tensor scales the model is unusable at 4 bits (\(\mathrm{PPL} = 3.03\), \(D = 1.60\) bits per token), marginal at 5 (\(\mathrm{PPL} = 1.068\)), and indistinguishable from fp64 at 6 (\(\mathrm{PPL} = 1.00023\), \(D = 2.1 \times 10^{-4}\) bits). Per-row scales buy that bit back: 4 bits already give \(\mathrm{PPL} = 1.114\) and 5 bits \(\mathrm{PPL} = 1.0027\). In rate terms, the eight matrices cost \(272 b\) bits — \(1360\) bits at the per-tensor knee (\(b = 5\)), \(1088\) at the per-row knee (\(b = 4\)) — against the \(133\) bits (\(32 \times \log_2 18\)) it takes to write the corpus down verbatim. The model does not pay for itself on 32 tokens at any bit width; the Information lens below turns that into a two-part code length.

Weights vs Activations

Theory

Weights are quantised once, offline, with a scale chosen from the whole tensor. Activations are produced at run time, token by token, and their statistics are worse in two ways. First, outlier channels: a few feature columns of the residual stream carry values many times larger than the rest, so a per-tensor scale sized for the outlier crushes the resolution of every other channel. Second, range varies per token, so a static scale is either too coarse for quiet tokens or saturates loud ones. The fixes mirror the weight case — a scale per token (per row of \(\mathbf{H}\)) chosen dynamically, or per channel with the outlier channels scaled into the weights (the SmoothQuant idea) — and the distortion that matters is again the KL of the softmax, now with the logits themselves on a \(b\)-bit grid.

Example — Quantising the residual stream and the logits

H = cache_fp.H_out;                      % residual stream entering the LM head, T x d_model
ch_max = zeros(d_model);
for j = 1:d_model
  ch_max(j) = max(abs(H(:, j)));
end
tok_max = max(abs(H), [], 2);
print("channel max|h|:", ch_max, "   per-token max|h| from", min(tok_max), "to", max(tok_max));
SNR_H = zeros(n_b, 2);
for i = 1:n_b
  [Hq, sq] = quant_tensor(H, bits(i));   SNR_H(i, 1) = snr_mat(H, Hq);
  [Hr, sr] = quant_rows(H, bits(i));     SNR_H(i, 2) = snr_mat(H, Hr);
end
print("H_out SNR at 8 bits: per-tensor", SNR_H(7, 1), " per-token", SNR_H(7, 2), "   at 4 bits:", SNR_H(3, 1), SNR_H(3, 2));
[lq4, s4] = quant_tensor(lo_fp, 4);                     % 4-bit logits
KL_logits4 = mean_kl(lo_fp, lq4, T - 1);
Pq = P;  Pq.W_U = quant_any(P.W_U, 4, 0);               % 4-bit W_U, everything else fp64
[lo_q, L_q, c] = transformer_forward(tokens, mask_train, Pq);
KL_WU4 = mean_kl(lo_fp, lo_q, T - 1);
print("KL at 4 bits — W_U only:", KL_WU4, "bits   logits:", KL_logits4, "bits   max|logit| =", max(max(abs(lo_fp))));
figure();
subplot(1, 2, 1)
bar(1:d_model, ch_max)
ylim([0, 1.1 * max(ch_max)])
title("Residual-stream channel ranges (H_out)")
xlabel("channel j"); ylabel("max |H(:, j)|")
subplot(1, 2, 2)
plot(bits, SNR_H(:, 1)', "color", "blue", "label", "per-tensor scale")
hold("on")
plot(bits, SNR_H(:, 2)', "color", "red", "label", "per-token scale")
hold("off")
title("SNR of the quantised residual stream")
xlabel("bits per activation"); ylabel("SNR (dB)")

channel max|h|: [1×4]  67.856716  26.318431  33.185321  24.311239    per-token max|h| from 10.66671248970273 to 67.85671611427112
H_out SNR at 8 bits: per-tensor 44.77927512063253  per-token 49.708883758554336    at 4 bits: 18.546316356983755 23.62533090710828
KL at 4 bits — W_U only: 0.0009081778263592612 bits   logits: 0.24446304315149164 bits   max|logit| = 117.07313909192767

The KL line is the point. Quantising the \(4 \times 18\) LM head to 4 bits costs \(\mathrm{KL} = 9.1 \times 10^{-4}\) bits per token; quantising the logits it produces to the same 4 bits costs \(0.24\) bits — a few hundred times more — because the logits span \(\pm 117\) with a handful of extreme entries, and a 4-bit step of that range is coarser than the gaps that decide the softmax. This is why "W4A16" (4-bit weights, 16-bit activations) is a common deployment point and "W4A4" is a research topic.

Saturate vs Wrap

Theory

An integer cast has to do something with a value outside its range. Saturation clips to the nearest rail; wraparound (two's-complement overflow) keeps the low bits and discards the carry, so 128 becomes \(-128\) and 200 becomes \(-56\). For a weight, saturation produces the closest representable value and a bounded error; wraparound flips the sign of the largest weights — precisely the ones that matter most — and the error is unbounded. Hardware wraps because it is free; a quantiser must saturate, and the demonstration below is the reason. The scenario is a stale calibration: a scale chosen when \(\lvert \mathbf{W}_U \rvert \le 0.3\) (the initialisation), applied after training stretched the range to \(1.58\).

Example — A stale scale under both overflow rules

W = P.W_U;
s_stale = 0.3 / 127;                                    % calibrated for |W| <= 0.3
q_sat  = int8(W / s_stale);                             % default: saturate
q_wrap = int8(W / s_stale, "wrap");                     % two's-complement overflow
print("overflowing entries:", sum(sum(abs(W / s_stale) > 127)), "of", numel(W));
print("W_U(1,1) =", W(1, 1), " -> W/s =", W(1, 1) / s_stale, "  saturate ->", q_sat(1, 1), "  wrap ->", q_wrap(1, 1));
snr_sat  = snr_mat(W, double(q_sat) * s_stale);
snr_wrap = snr_mat(W, double(q_wrap) * s_stale);
Pq = P;  Pq.W_U = double(q_sat) * s_stale;
[lo_q, L_sat, c] = transformer_forward(tokens, mask_train, Pq);
Pq = P;  Pq.W_U = double(q_wrap) * s_stale;
[lo_q, L_wrap, c] = transformer_forward(tokens, mask_train, Pq);
print("saturate: SNR", snr_sat, "dB  PPL", exp(L_sat), "     wrap: SNR", snr_wrap, "dB  PPL", exp(L_wrap));

overflowing entries: 59 of 72
W_U(1,1) = 0.4330597614525455  -> W/s = 183.32863234824427   saturate -> 127   wrap -> -73
saturate: SNR 4.04437859575337 dB  PPL 4.59455301907125      wrap: SNR -0.6855708294558249 dB  PPL 967932610.3346702

\(59\) of the 72 entries overflow. Saturating them costs a lot — SNR \(4.0\) dB, \(\mathrm{PPL} = 4.59\) — but the model still ranks tokens sensibly. Wrapping turns \(W_U(1,1) = 0.433\) into the integer \(-73\), a sign flip; the SNR goes negative (\(-0.7\) dB: the error is larger than the signal) and perplexity explodes to \(\mathrm{PPL} \approx 9.7 \times 10^{8}\). Overflow should never be silent.

KV-Cache Memory

Theory

During generation the KV cache (21-sampling-and-generation) stores one key and one value vector per layer, per KV head, per past token. With \(N_{\text{layers}}\) layers, \(H_{\text{kv}}\) key/value heads of width \(d_{\text{head}}\), and \(b\) bits per element,

\[\text{bytes per token} = 2 \cdot N_{\text{layers}} \cdot H_{\text{kv}} \cdot d_{\text{head}} \cdot \frac{b}{8}, \qquad \text{cache} = T \times \text{bytes per token}.\]

GQA (24-modern-architectural-variants) shrinks \(H_{\text{kv}}\); quantisation shrinks \(b\); the two multiply. The cache is also read once per decoded token, so it joins the weights in the denominator of the bandwidth bound — Exact as an upper bound on tokens per second.

Example — LLaMA-2-70B at 16, 8 and 4 bits

N_layers = 80;  H_kv = 8;  d_head = 128;                % LLaMA-2-70B (GQA)
T_ctx = 8192;
kv_bits = [16, 8, 4];
bytes_tok = 2 * N_layers * H_kv * d_head * kv_bits / 8;
cache_GB = bytes_tok * T_ctx / 1e9;
tok_both = BW ./ ((w_GB + cache_GB) * 1e9);             % weights and cache at the same width
print("bits:                 ", kv_bits);
print("KiB per token:        ", bytes_tok / 1024);
print("cache at T = 8192 (GB):", cache_GB);
print("weights (GB):         ", w_GB);
print("tok/s bound (both):   ", tok_both);
T_sweep = (1:32) * 4096;
figure();
subplot(1, 2, 1)
hold("on")
plot(T_sweep / 1024, bytes_tok(1) * T_sweep / 1e9, "color", "blue", "label", "16-bit cache")
plot(T_sweep / 1024, bytes_tok(2) * T_sweep / 1e9, "color", "red", "label", "8-bit cache")
plot(T_sweep / 1024, bytes_tok(3) * T_sweep / 1e9, "color", "green", "label", "4-bit cache")
yline(80, "gray", "one 80 GB device")
hold("off")
title("LLaMA-2-70B KV cache vs context length")
xlabel("context length (k tokens)"); ylabel("GB")
subplot(1, 2, 2)
bar({"16-bit", "8-bit", "4-bit"}, tok_both)
legend("tok/s upper bound")
ylim([0, 100])
title("Bandwidth-bound decode rate (weights + cache)")
xlabel("bits per weight and per cache element"); ylabel("tok/s upper bound")

bits:                  [1×3]  16.000000  8.000000  4.000000
KiB per token:         [1×3]  320.000000  160.000000  80.000000
cache at T = 8192 (GB): [1×3]  2.684355  1.342177  0.671089
weights (GB):          [1×3]  140.000000  70.000000  35.000000
tok/s bound (both):    [1×3]  23.478398  46.956795  93.913590

Engineering Lenses

Signals

Exact. Round-to-nearest quantisation is an additive error, \(\hat{\mathbf{W}} = \mathbf{W} + \mathbf{e}\) with \(\lvert\mathbf{e}\rvert \le \Delta/2\), and given an error uniform on \(\pm\Delta/2\) its variance is \(\Delta^2/12\). Model. Treating that error as signal-independent uniform noise — the uniform-noise model of the ADC — is what turns \(\Delta^2/12\) into 6.02 dB of SNR per bit; the weight tensor plays the role of the sampled waveform and \(\max\lvert\mathbf{W}\rvert\) that of the full-scale range, and the SNR sweep above already showed the model breaking at 2–3 bits. The pooled errors of all 272 quantised weights, normalised by their steps, test the model directly. The second computation is the rounding-mode corollary: qfmt's default mode is floor (truncation), whose error is uniform on \([0, \Delta)\) with mean square \(\Delta^2/3\) — exactly 6.02 dB worse than rounding to nearest. It uses \(\mathcal{N}(0, 0.3^2)\) samples in Q1.6, the format that holds their range; in Q0.7 the few samples beyond \(\pm 1\) clip, and clipping error — outside the uniform-noise model — then dominates the SNR (which is why quoted Q0.7 figures for this signal vary by several dB with the seed).

b_sig = 4;
e_pool = [];
for k = 1:8
  W = mat_k(P, k);
  [Wt, st] = quant_tensor(W, b_sig);
  e_pool = [e_pool, reshape(W - Wt, 1, numel(W)) / st];     % error in units of the step
end
print("4-bit errors e/Delta over", length(e_pool), "weights: mean", mean(e_pool), "  12 var", 12 * mean(e_pool .^ 2), "  max|e|/Delta", max(abs(e_pool)));
seed(26);
x = randn(10000) * 0.3;                                   % Gaussian "weights", max |x| ~ 1.4 < 2
snr_round = snr(x, quantize(x, qfmt(8, 6, "round", "saturate")));
snr_floor = snr(x, quantize(x, qfmt(8, 6)));              % default: floor rounding
snr_model = 10 * log10(0.09 / ((2 ^ -6) ^ 2 / 12));
print("N(0, 0.3^2) in Q1.6:  round-to-nearest", snr_round, "dB (model", snr_model, ")   truncate", snr_floor, "dB   difference", snr_round - snr_floor, "dB");
figure();
histogram(e_pool, 12);
title("Normalised quantisation error e / Delta of the 272 weights at 4 bits")
xlabel("e / Delta"); ylabel("count")

4-bit errors e/Delta over 272 weights: mean -0.0134772696577325   12 var 0.9721392833809697   max|e|/Delta 0.4971076506487199
N(0, 0.3^2) in Q1.6:  round-to-nearest 36.42951797969648 dB (model 36.45783703454724 )   truncate 30.392707555327853 dB   difference 6.036810424368625 dB

Systems

Model. Post-training quantisation, as done above, is open-loop: the quantiser is applied after the training loop closed. Quantisation-aware training puts the quantiser inside the loop — the forward pass sees \(\hat{\mathbf{W}} = s\,\operatorname{round}(\mathbf{W}/s)\), the backward pass treats the rounding as the identity (the straight-through estimator, since \(\partial \operatorname{round}/\partial W\) is zero almost everywhere), and the optimiser updates a full-precision master copy of \(\mathbf{W}\). The loop then steers the master weights toward values whose rounded versions fit the corpus, instead of hoping the rounded versions of the fp optimum do. A few lines turn the capstone trainer into that loop; it starts from the 3-bit per-tensor model, whose post-training perplexity was \(14.7\).

b_qat = 3;
Pm = P;                                                  % full-precision master copy
[Mq, Vq] = adamw_init(Pm);
Pq = quant_params(Pm, b_qat, 0);
[lo_q, L_before, c] = transformer_forward(tokens, mask_train, Pq);
for k = 1:100
  Pq = quant_params(Pm, b_qat, 0);                       % quantiser inside the forward pass
  [lo_q, L_q, cache_q] = transformer_forward(tokens, mask_train, Pq);
  dl = ce_dlogits(tokens, mask_train, lo_q, cache_q.total);
  G = transformer_backward(tokens, dl, cache_q, Pq);     % gradient at the quantised weights ...
  [Pm, Mq, Vq] = adamw_step(Pm, G, Mq, Vq, 0.002, k, bet1, bet2, eps_a, 0.0);   % ... applied to the master (STE)
end
Pq = quant_params(Pm, b_qat, 0);
[lo_q, L_after, c] = transformer_forward(tokens, mask_train, Pq);
print("3-bit per-tensor:  post-training PPL", exp(L_before), "  ->  after 100 QAT steps", exp(L_after), "   KL", mean_kl(lo_fp, lo_q, T), "bits");

3-bit per-tensor:  post-training PPL 14.725295329311486   ->  after 100 QAT steps 1.0112351775253505    KL 0.015316221039036801 bits

The 3-bit model goes from \(\mathrm{PPL} = 14.7\) to \(1.011\) without changing the bit width — the rate stayed at \(272 \times 3\) bits and the loop found a better codebook entry. The label is Model rather than Exact because the straight-through estimator is not the gradient of the quantised loss (that gradient is zero); it is a surrogate that works in practice.

Information

Exact. A quantised model is a code for the corpus in the two-part sense of MDL (23-putting-it-all-together): first transmit the weights (\(272\,b\) bits), then the corpus under the model (\((T - 1)\,\mathcal{L}_q / \ln 2\) bits — \(\mathcal{L}_q\) is a mean over the \(T - 1 = 31\) predicted tokens). The rate–distortion curve measured above is the trade between the two parts, and their sum is the description length to minimise.

bits_corpus_raw = T * log2(vocab);
code_model = n_q * bits;
code_corpus = (T - 1) * L_t / log(2);
two_part = code_model + code_corpus;
for i = 1:n_b
  print(sprintf("%d bits/weight:  weights %5d + corpus %8.2f = %8.1f bits", bits(i), code_model(i), code_corpus(i), two_part(i)));
end
[best_len, best_i] = min(two_part);
print("raw corpus:", bits_corpus_raw, "bits   fp64 model:", n_q * 64, "bits   best two-part code:", best_len, "bits at b =", bits(best_i));

2 bits/weight:  weights   544 + corpus   416.90 =    960.9 bits
3 bits/weight:  weights   816 + corpus   120.29 =    936.3 bits
4 bits/weight:  weights  1088 + corpus    49.62 =   1137.6 bits
5 bits/weight:  weights  1360 + corpus     2.96 =   1363.0 bits
6 bits/weight:  weights  1632 + corpus     0.01 =   1632.0 bits
7 bits/weight:  weights  1904 + corpus     0.01 =   1904.0 bits
8 bits/weight:  weights  2176 + corpus     0.00 =   2176.0 bits
raw corpus: 133.437600046154 bits   fp64 model: 17408 bits   best two-part code: 936.2869645779949 bits at b = 3

The description length is minimised at \(b = 3\) (\(936\) bits) — below the perplexity knee, because once the corpus term is small the weight term dominates — and every choice of \(b\) exceeds the \(133\) bits of writing the 32 tokens down raw. That is the honest information-theoretic reading of the capstone: 300 parameters memorised 133 bits of corpus. A real language model earns its bits because the corpus term is \(10^{13}\) tokens long and the weight term is not; the shape of the trade — a knee in the corpus term, a linear weight term — is the same.

Key Takeaways

• Symmetric uniform quantisation is \(s = \max\lvert\mathbf{W}\rvert / (2^{b-1}-1)\), \(\mathbf{q} = \operatorname{clip}(\operatorname{round}(\mathbf{W}/s))\), \(\hat{\mathbf{W}} = s\mathbf{q}\); a per-tensor scale is a Q-format whose step is tuned to the data, and per-row scales buy up to one bit (8 dB on \(\mathbf{E}\), 3.5 dB averaged over the eight matrices).
• Quantisation error is additive uniform noise with variance \(\Delta^2/12\): 6.02 dB per bit minus a crest-factor penalty, measured here at \(6.31\) dB per bit from 5 to 8 bits, and breaking down below 4 bits.
• The rate–distortion curve of the capstone model has its knee at 5–6 bits per weight with per-tensor scales and 4–5 with per-row scales; KL of the next-token distribution is the distortion to report.
• Activations are harder than weights: outlier channels and per-token range make 4-bit logits cost hundreds of times the KL of 4-bit weights.
• Saturate, never wrap: a stale scale under wraparound turns the largest weights into their negatives and destroys the model.
• KV-cache bytes per token are \(2 N_{\text{layers}} H_{\text{kv}} d_{\text{head}}\, b / 8\); decode rate is bounded by bandwidth over bytes read per token, and halving \(b\) for weights and cache together doubles it.

Standalone Scripts

Run with make lesson-26 or rustlab run lessons/26-quantization-and-fixed-point-inference/<script>.rlab.

Expected Numerical Outputs Summary

Exercises

1. Per-column scales. quant_rows puts one scale per input channel of a projection. Write quant_cols with one scale per output column and show algebraically that \(\mathbf{X}\hat{\mathbf{W}} = (\mathbf{X}\mathbf{q}) \operatorname{diag}(\mathbf{s})\) — the scales factor out of the integer accumulation, which is why deployed kernels prefer that axis. Compare SNR and PPL at 4 bits against the per-row and per-tensor curves.

2. Derive \(\Delta^2/12\) and the 6.02 dB. Starting from an error uniformly distributed on \([-\Delta/2, \Delta/2]\), derive its variance, then show that adding one bit changes the SNR by \(20\log_{10} 2\) dB. Where does the \(+4.77\) dB in the SNR formula come from, and why does a Gaussian tensor with \(\max\lvert W\rvert = 4\sigma\) land at \(6.02\,b - 7.3\) dB?

   Solution

   For \(e \sim \mathcal{U}[-\Delta/2, \Delta/2]\) the mean is zero by symmetry and the variance is \(\sigma_e^2 = \frac{1}{\Delta}\int_{-\Delta/2}^{\Delta/2} e^2\,de = \frac{1}{\Delta}\left[\frac{e^3}{3}\right]_{-\Delta/2}^{\Delta/2} = \frac{1}{\Delta}\cdot\frac{2}{3}\cdot\frac{\Delta^3}{8} = \frac{\Delta^2}{12}\). With full-scale range \(R = \max\lvert W\rvert\) and \(2^{b-1}\) codes per side, \(\Delta \approx R\,2^{1-b}\), so \(\sigma_e^2 = R^2\,2^{2-2b}/12 = R^2\,2^{-2b}/3\). Adding a bit multiplies \(\sigma_e^2\) by \(1/4\): \(10\log_{10}4 = 20\log_{10}2 = 6.02\) dB. The SNR is \(10\log_{10}\!\big(P_W \cdot 3 \cdot 2^{2b}/R^2\big) = 6.02\,b + 10\log_{10}3 - 20\log_{10}(R/\sqrt{P_W}) = 6.02\,b + 4.77 - 20\log_{10}c\); the \(4.77\) is the \(\tfrac{1}{3}\) in \(\sigma_e^2\). For a Gaussian tensor clipped at \(R = 4\sigma\), \(c = 4\) and \(20\log_{10}4 = 12.04\), giving \(6.02\,b - 7.3\) dB. (The textbook ADC figure \(6.02\,b + 1.76\) dB is the same formula for a full-scale sine, whose crest factor is \(\sqrt{2}\): \(4.77 - 3.01 = 1.76\).)

3. Mixed precision. Quantise \(\mathbf{E}\) and \(\mathbf{W}_U\) at 8 bits and the six inner matrices at 4 bits (or the reverse) and measure PPL and KL. Which matrices are most sensitive, and does the answer agree with their crest factors? The "only \(\mathbf{E}\)" and "only \(\mathbf{W}_U\)" sweeps are two-line changes to quant_params.

4. Smooth the outlier channel. Divide column 1 of H by \(\alpha = $ `ch_max(1) / ch_max(2)` and multiply row 1 of $\mathbf{W}_U\) by the same \(\alpha\) — the product \(\mathbf{H}\mathbf{W}_U\) is unchanged. Re-measure the 4-bit SNR of the smoothed \(\mathbf{H}\) and the KL of the 4-bit logits. This is the SmoothQuant transformation.

5. Where the cache overtakes the weights. From bytes_tok, solve for the context length at which the KV cache equals the weight bytes, at fp16 for both and at 4 bits for both. Then redo the decode bound at 128k tokens: which term dominates, and what does GQA with \(H_{\text{kv}} = 8\) versus MHA with \(H_{\text{kv}} = 64\) do to it?

What's next

This is the last lesson. The course began in 00-the-llm-as-a-system with a picture — a chain of operators on a \(T \times d\) signal inside two feedback loops — and three questions to ask of every box; you have now asked them of every box, from the alphabet to the bit width of the deployed weights. The real systems are within reach: nanoGPT's model.py is 14-full-gpt-architecture in PyTorch; llama.cpp's k-quants are the per-block scales of this lesson with 2–6 bits and a second-level scale for the scales; GPTQ and AWQ are post-training quantisers that choose the rounding and the per-channel scales to minimise the output error rather than the weight error — the KL panel of this lesson turned into an objective. Every one of them is built from the pieces in this course.