Model·Foundations

What the Hessian means in a quantization paper

Quantization and pruning papers lean on 'the Hessian' constantly, and they don't all mean the same matrix. Notes on what it is, why the layer-wise version is exactly computable, and the three different jobs it does across the literature. With small numerical experiments, one of which I got wrong on the first try.

Note
concept reviewed published July 12, 2026 · updated July 18, 2026

Core insights. The Hessian is the matrix of second derivatives of a loss with respect to the weights. For a network’s full training loss it is astronomically large and everyone approximates it. But for the layer-wise reconstruction objective that post-training quantization actually optimizes, the Hessian is small, exact, and cheap: H=2XXH = 2XX^\top, where XX is the calibration input. Its diagonal says which weights are dangerous to round, and its off-diagonal says how to fix the damage using other weights. Different methods use those two pieces of information differently, and once I sorted the literature by which Hessian and which job, papers that looked similar stopped looking similar.

Second derivatives, in one paragraph

Take a loss ff over weights wRdw \in \mathbb{R}^d. The gradient tells you the slope; the Hessian Hij=2f/wiwjH_{ij} = \partial^2 f / \partial w_i \partial w_j tells you the curvature: how quickly the loss changes if you move weights ii and jj together. Around a point where the gradient is zero, Taylor expansion gives

f(w+δ)12δHδ,f(w + \delta) \approx \tfrac{1}{2}\,\delta^\top H\, \delta ,

so the Hessian is precisely the object that prices a perturbation. Quantization is a perturbation you’re forced to make. That is the whole reason this matrix keeps appearing: every “how much does rounding hurt and where” question is a question about HH.

The layer-wise Hessian, derived

Papers about second-order methods on the training loss have to fight for their Hessian: for a 7B model it’s a 7B×7B7\text{B} \times 7\text{B} matrix that nobody can even store. The move that makes post-training quantization tractable is changing the question. Don’t preserve the training loss; preserve each layer’s outputs on calibration data.

For one output neuron (one row ww of the layer’s weight matrix) and calibration inputs XRd×nX \in \mathbb{R}^{d \times n}:

f(w^)=wXw^X22=δ(XX)δ,δ=w^w.f(\hat w) = \lVert wX - \hat wX \rVert_2^2 = \delta^\top (X X^\top)\, \delta , \qquad \delta = \hat w - w .

Two derivatives give H=2XXH = 2XX^\top. Note what happened: this objective is not approximately quadratic near a minimum, it is exactly quadratic everywhere, so the Taylor expansion has no error term. And HH is d×dd \times d where dd is the layer’s input width, a few thousand, not a few billion.

I checked this numerically rather than trusting my own algebra: finite-difference second derivatives of ff on random data (d=8d=8, 200 samples) match 2XX2XX^\top to a relative error of 1.2×1071.2 \times 10^{-7}, which is finite-difference noise.

Two structural facts fall out of the formula. First, XX is the only ingredient, so every row of the layer shares the same Hessian; each neuron’s error surface is the same paraboloid translated to a different center. GPTQ’s speed comes from this: one Cholesky per layer, reused across all rows. Second, HH depends on the calibration data, not on the weights. Your Hessian is only as representative as your 128 samples of C4.

Reading the matrix

The diagonal entry is Hii=2txi(t)2H_{ii} = 2\sum_t x_i(t)^2, the energy of input feature ii over the calibration set. A weight multiplying a loud, frequent input is expensive to round; a weight on a quiet input is nearly free. When GPTQ’s actorder sorts columns by diag(H) descending, this is the quantity it sorts by.

The off-diagonal entry is Hij=2txi(t)xj(t)H_{ij} = 2\sum_t x_i(t)\,x_j(t), the correlation between inputs ii and jj, and it is what makes error compensation possible rather than just error avoidance. The degenerate case makes it obvious: if x2x_2 is always equal to x1x_1, then y=w1x1+w2x2=(w1+w2)x1y = w_1 x_1 + w_2 x_2 = (w_1 + w_2)x_1, and any rounding error on w1w_1 can be exactly cancelled by shifting w2w_2. I ran exactly this: forcing an error of 0.3 onto w1w_1 produced a reconstruction error of 16.6, and subtracting 0.3 from w2w_2 brought it to 102910^{-29}. Zero, in floats.

For general correlations, Optimal Brain Surgeon gives the optimal version of that trick. Freeze wqw_q at its rounded value, and the best adjustment to the remaining weights FF is

δF=wqquant(wq)[HF1]qq  (HF1):,q.\delta_F = -\,\frac{w_q - \mathrm{quant}(w_q)}{[H_F^{-1}]_{qq}}\;(H_F^{-1})_{:,q} .

I verified this against brute force: quantize one coordinate, apply the OBS update, and separately solve the constrained least-squares problem exactly. The two weight vectors agree to 7×10167 \times 10^{-16} and the losses to six digits. An embarrassing detail that belongs in this note: my first attempt flipped the sign of the error term, and “compensation” made the loss worse than doing nothing. The brute-force cross-check is what caught it. If you implement any of these formulas, build the dumb exact baseline first; sign errors in compensation code fail silently and look like “the method just doesn’t help much.”

With the sign right, I ran an end-to-end miniature: 8 correlated inputs, a coarse 0.25 grid, sequential quantize-and-compensate against round-to-nearest, over 500 random instances. Median result: compensation is 1.46× better, winning about 70% of instances, with a long tail of big wins (up to 22×) and real losses on unlucky draws (down to 0.21×). My first draft quoted a single 5.3× instance here; when a reshuffled random stream later produced a losing draw, I realized the single number was a cherry-pick and replaced it with the distribution. The lesson is worth stating: the greedy sequence is locally optimal at every step but carries no per-instance guarantee at toy scale, and the evidence that it wins decisively at real scale is the paper’s own OPT and BLOOM tables, not my miniature.

The same matrix, three different jobs

Sorting the literature by what the Hessian is for:

Job 1: decide what is sensitive. Optimal Brain Damage (LeCun, 1990) started this, pruning weights by saliency 12wi2Hii\tfrac{1}{2} w_i^2 H_{ii} with a diagonal approximation of the training-loss Hessian. HAWQ and HAWQ-V2 (Dong, Yao et al., 2019) are the modern descendants: they never compensate anything, they use the loss Hessian’s top eigenvalue (V1) or its trace, estimated by Hutchinson sampling (V2), to rank whole layers by fragility and assign bit-widths, high curvature getting more bits. SqueezeLLM (2023) does the per-weight version for LLMs, using a diagonal Fisher-information approximation to pull k-means centroids toward sensitive weights.

Job 2: fix the damage. Optimal Brain Surgeon (Hassibi & Stork, 1992) added the inverse-Hessian compensation update above, arguing OBD’s diagonal assumption throws away the useful part. Optimal Brain Compression (Frantar & Alistarh, 2022) made OBS practical per layer with the exact 2XX2XX^\top Hessian, and GPTQ is that machinery made fast enough for 175B models. The compensation family needs off-diagonals by definition; a diagonal Hessian would leave it nothing to work with.

Job 3: decide which way to round. AdaRound (Nagel et al., 2020) is the one I think people misfile. Its Taylor analysis shows that round-to-nearest is only optimal if the Hessian were diagonal; with off-diagonals, rounding some weights up and some down in a coordinated pattern beats nearest rounding. It then optimizes the binary up/down choices against the same layer-wise E[xx]E[xx^\top] objective, weights staying on grid points rather than moving freely. So GPTQ and AdaRound optimize the same quadratic with different degrees of freedom: continuous compensation versus binary rounding direction.

The split that matters when reading any of these papers: which loss. OBD, HAWQ, and Fisher methods differentiate the end-to-end training loss (huge, approximated, needs labels or at least loss evaluations). OBC, GPTQ, and AdaRound in practice differentiate the layer reconstruction loss (small, exact, needs only activations). When a paper says “Hessian-aware,” check which one before assuming anything transfers.

BRECQ and the unit of reconstruction

That split between the two losses has a middle I left out, and it is where BRECQ (Li et al., 2021) lives. Layer reconstruction preserves one layer’s outputs at a time; end-to-end methods preserve the whole training loss. BRECQ reconstructs a block: several consecutive layers at once, a residual bottleneck in its own ResNet and MobileNet experiments, or a transformer block once the idea is carried to LLMs. The reason to want a unit this size is the worry this note keeps circling back to: a per-layer objective is blind to everything downstream of its own layer, and a whole-network objective overfits a calibration set of a thousand images. A block is the compromise.

The math that licenses it is worth walking through slowly, because it is the same move that made the layer-wise Hessian legal, applied one level up. What you actually want to minimize is the increase in task loss when rounding perturbs the weights by Δw\Delta w. Expand it to second order. Because the pretrained model already sits at a minimum in weight space, the gradient term of BRECQ’s Eq. (4) is close to zero, and the second-order piece is what remains:

E[L(w+Δw)]E[L(w)]12ΔwH(w)Δw.\mathbb{E}[L(w + \Delta w)] - \mathbb{E}[L(w)] \approx \tfrac{1}{2}\,\Delta w^\top H^{(w)} \Delta w .

H(w)H^{(w)} is the task-loss Hessian for a whole block’s worth of weights, far too large to form. The Gauss-Newton approximation (Eq. 6) rewrites any such weight-Hessian through the Jacobian of an output: with J=z/wJ = \partial z / \partial w, H(w)JH(z)JH^{(w)} \approx J^\top H^{(z)} J. Apply that to a block’s own output zz and substitute, and the quadratic you cannot build collapses into one you can measure:

ΔwH(w)Δw(JΔw)H(z)(JΔw)=ΔzH(z)Δz.\Delta w^\top H^{(w)} \Delta w \approx (J\Delta w)^\top H^{(z)} (J\Delta w) = \Delta z^\top H^{(z)} \Delta z .

Minimizing a weight-space object you cannot store became minimizing the block’s output error Δz\Delta z, weighted by the output Hessian H(z)H^{(z)}. One honest gap in that last step: JΔwJ\Delta w is the linear part of the output change, and BRECQ actually minimizes the true quantized-minus-original block output, which equals JΔwJ\Delta w only to first order once the block has nonlinearities in it. This is the same trick behind the H=2XXH = 2XX^\top from earlier. For a single linear layer with z=wXz = wX the Jacobian is XX^\top, and plain output MSE is the case H(z)=2IH^{(z)} = 2I, so JH(z)J=2XXJ^\top H^{(z)} J = 2XX^\top exactly. I checked that equality, and checked that a non-identity H(z)H^{(z)} gives a genuinely different weight-space geometry: a Fisher reweighting whose entries span orders of magnitude leaves a 28% residual against the best rescaling of 2XX2XX^\top. So the layer-wise Hessian this whole note is built on is the block-reconstruction Hessian with the block shrunk to one layer and the output metric set to the identity.

That output metric is BRECQ’s other idea. H(z)H^{(z)}, the task loss curved against the block’s own output, is itself dense and unknown, so BRECQ approximates it by its diagonal, and approximates that by the empirical Fisher: the squared task gradients diag((L/zi)2)\operatorname{diag}\big((\partial L/\partial z_i)^2\big). The intuition is concrete. An output coordinate the final loss is sensitive to (large L/zi\partial L / \partial z_i) should have its reconstruction error punished harder than a coordinate the loss barely notices. So the per-block objective (Eq. 10) is igi2Δzi2\sum_i g_i^2\,\Delta z_i^2 averaged over calibration data, with gi=L/zig_i = \partial L/\partial z_i. Layer reconstruction, in GPTQ and AdaRound, weights every output coordinate equally; BRECQ weights them by how much the task cares. I checked the identity underneath this on a softmax head, where the Gauss-Newton Fisher equals the output Hessian exactly, H(z)=diag(p)ppH^{(z)} = \operatorname{diag}(p) - pp^\top. There the squared-gradient diagonal is an unbiased estimate of diag(H(z))=pi(1pi)\operatorname{diag}(H^{(z)}) = p_i(1-p_i) when the label is drawn from the model, but BRECQ uses the data’s label and each calibration example gives one noisy draw (up to 3× noisier than its own mean in my toy). That noise is the empirical-Fisher asterisk from the next section, showing up as the actual object BRECQ optimizes.

Then I read the code, and the theory and the default part company. The Fisher-weighted objective is implemented (fisher_diag, block_recon.py:150-151, commit e455d62) as literally ((pred - tgt).pow(2) * grad.pow(2)).sum(1).mean(), which is Eq. (10) transcribed. But the released entry point never turns it on: main_imagenet.py:202 passes opt_mode='mse', so the reconstructions that reproduce the paper’s tables minimize unweighted output error (H(z)=IH^{(z)} = I), and the gradients Eq. (10) needs are never even cached (block_recon.py:78). The Fisher apparatus the paper derives is off by default in the code that runs it. There is also a third mode, fisher_full, that approximates the off-diagonal Fisher with a per-sample outer product divided by a bare 100 (block_recon.py:152-156), which I do not find in the paper’s equations. I read the code and did not re-run ImageNet, so I cannot say how much the weighting moves accuracy; I am reporting the default, not claiming MSE is worse.

The last choice, why a block and not a layer or the network, BRECQ argues as bias against variance: the block-diagonal Hessian keeps the cross-layer coupling inside a block and discards the coupling between blocks, which lowers approximation error while limiting how much the fit can overfit the calibration set. I did not verify this myself; it needs a real network for the overfitting to appear. The cleanest evidence I have is second-hand and at a different scale, in the EfficientQAT note, where reconstructing one block at a time beats both per-layer and multi-block units on a 2-bit LLM (5.93, against 6.10 per-layer and 6.00 for four blocks). So “block-diagonal” is not a claim that the Hessian is block-diagonal. It is a regularization knob, and where the knob sits best is measured, not derived.

Common confusions

The Fisher-equals-Hessian shortcut deserves its asterisk. The Fisher information matrix equals the expected Hessian under the model’s own predictive distribution; the empirical Fisher that methods like SqueezeLLM actually compute (squared gradients on data) is a further approximation, and it degrades when the model is far from fit or gradients are near zero. It’s a reasonable sensitivity proxy and a poor curvature estimate, which is fine for Job 1 and would be a disaster for Job 2.

Damping is not a detail. XXXX^\top is only positive semi-definite: any input direction the calibration set never excites has zero curvature, and the inverse doesn’t exist. Every real implementation adds λI\lambda I before inverting (GPTQ uses 1% of the mean diagonal, gptq.py:98–100). Why one percent is reliably enough is a question I have not seen analyzed anywhere; every repo just inherits the constant.

And the one I’d tattoo somewhere visible: the layer-wise Hessian knows nothing about the network past its own layer. “Preserve every layer’s outputs” is a proxy for “preserve the model,” and the gap between those two objectives is exactly where layer-local methods lose capability that perplexity doesn’t show. BRECQ’s block unit, above, is the most direct attempt I know of to narrow that gap without paying for full end-to-end optimization, and it only narrows it: a block still cannot see past its own last layer.

What I verified

The numerical checks in this note (finite-difference Hessian, exact compensation under perfect correlation, OBS-equals-brute-force, the ρ closed form, the 500-instance RTN-versus-compensation experiment, the Gauss-Newton identity that makes 2XX2XX^\top a special case of BRECQ’s block Hessian, and the softmax-head fact that the empirical Fisher diagonal is a one-sample estimate of the output-Hessian diagonal) are one seeded NumPy script, published in the site’s repo at verification/hessian-in-quantization/checks.py; run it and the numbers quoted here print as PASS lines. The GPTQ code references are from my audit of the reference repo, and the BRECQ references are from reading its reference repo at commit e455d62; the opt_mode='mse' default and the fisher_diag objective are both quoted from that source. BRECQ’s block-granularity argument I did not verify here; that one block wins is measured in the EfficientQAT note, not in this script. The characterizations of OBD, OBS, HAWQ, AdaRound, and SqueezeLLM are from reading those papers’ method sections, not their code; if any of their implementations diverge from the papers the way GPTQ’s does, I wouldn’t know yet.