Iterative Hessian Mixing
12 min read
Publishing the coefficients of a least squares fit tells the world something about every row that went into it, and differential privacy asks for the fit anyway, with enough noise that no one row is detectable. Lev, Shenfeld, Srinivasan, Ligett & Wilson [LSSLW, ICML ’26] answer with an algorithm that was already private before anyone asked it to be.
They take the iterative Hessian sketch of Pilanci & Wainwright [PW, JMLR ’16], a randomized solver built for speed a decade before anyone thought to privatize it, and observe that its two moving parts are, almost unchanged, two differential privacy mechanisms. The random sketch already releases Gaussian noise. The exact gradient step is a bounded sum, ready for the classical Gaussian mechanism. The resulting algorithm, iterative Hessian mixing (IHM), is the current state of the art for private least squares.
The Problem
The dataset is a matrix \(\mathbf{X} \in \mathbb{R}^{n \times d}\) of covariate rows \(\mathbf{x}_i \in \mathbb{R}^d\) and a response vector \(\mathbf{y} \in \mathbb{R}^n\), one row per person, normalized so that \(\|\mathbf{x}_i\|_2 \le 1\) and \(|y_i| \le 1\). Least squares seeks \(\boldsymbol{\theta}^* = \arg\min_{\boldsymbol{\theta}} L(\boldsymbol{\theta})\) for \(L(\boldsymbol{\theta}) = \|\mathbf{y} - \mathbf{X}\boldsymbol{\theta}\|_2^2\), and an algorithm’s error is its excess risk \(L(\widehat{\boldsymbol{\theta}}) - L(\boldsymbol{\theta}^*)\). The privacy constraint is \((\varepsilon, \delta)\)-differential privacy against replacing one person’s row by zeros: blanking out that row moves the probability of any output event by at most a factor \(e^{\varepsilon}\) plus a slack \(\delta\). (The standard tool is Gaussian noise scaled to the query’s sensitivity; we lean throughout on its privacy loss.)
The standard baseline is AdaSSP, due to Wang [Wang, UAI ’18]: add Gaussian noise to the sufficient statistics \(\mathbf{X}^\top\mathbf{X}\) and \(\mathbf{X}^\top\mathbf{y}\), add an adaptive ridge term sized to a privately released smallest eigenvalue, and solve. It is simple and matches the information-theoretic lower bounds up to logarithmic factors. But its risk bound carries a multiplier we will call \(\gamma_{\mathrm{a}}\):
\[ \gamma_{\mathrm{a}} = \frac{\sqrt{d \,\ln(1/\varrho) \,\ln(1/\delta)}}{\varepsilon}, \]
where \(\varrho\) is the probability that the bound is allowed to fail. When \(\varrho\) is kept exponentially small, that multiplier can sit a factor \(\sqrt{d}\) above the lower bound. Removing it is what IHM is for.
The Engine
IHM privatizes a solver that already existed, so we start with the solver. We expand the objective:
\[ L(\boldsymbol{\theta}) = \|\mathbf{y}\|_2^2 + \|\mathbf{X}\boldsymbol{\theta}\|_2^2 - 2\,\boldsymbol{\theta}^\top\mathbf{X}^\top\mathbf{y}. \]
The data enters through a quadratic term, whose curvature is the Gram matrix \(\mathbf{X}^\top\mathbf{X}\), and a linear term, the gradient’s data half. The classic sketch-and-solve recipe compresses everything at once: draw \(\mathbf{S} \in \mathbb{R}^{k \times n}\) with i.i.d. standard Gaussian entries and \(k \ll n\), then solve the \(k\)-row problem on \((\mathbf{S}\mathbf{X}, \mathbf{S}\mathbf{y})\). Its excess risk scales with the residual \(L(\boldsymbol{\theta}^*)\), discounted only by roughly \(\sqrt{d/k}\), so driving the error down requires a larger, less compressed sketch. Pilanci & Wainwright’s iterative Hessian sketch (IHS) instead sketches only the quadratic term, keeps the gradient exact, and iterates with a fresh, independent sketch \(\mathbf{S}_t\) each round, so a bad draw costs one round rather than the run:
\[ \widehat{\boldsymbol{\theta}}_{t+1} = \widehat{\boldsymbol{\theta}}_t + \Big(\tfrac{1}{k}(\mathbf{S}_t\mathbf{X})^\top(\mathbf{S}_t\mathbf{X})\Big)^{-1} \mathbf{X}^\top\big(\mathbf{y} - \mathbf{X}\widehat{\boldsymbol{\theta}}_t\big). \]
This is Newton’s method with the true gradient and a cheap, crude Hessian. Because \(\mathbf{S}_t\) has i.i.d. standard Gaussian entries, its scaled Gram matrix averages to the identity:
\[ \mathbb{E}\Big[\tfrac{1}{k}\mathbf{S}_t^\top\mathbf{S}_t\Big] = \mathbf{I}_n. \]
Sandwiching \(\mathbf{X}\) on both sides makes the sketched curvature unbiased for the Gram matrix:
\[ \mathbb{E}\Big[\tfrac{1}{k}(\mathbf{S}_t\mathbf{X})^\top(\mathbf{S}_t\mathbf{X})\Big] = \mathbf{X}^\top\mathbf{X}. \]
A crude Hessian costs speed and nothing else: at \(\widehat{\boldsymbol{\theta}}_t = \boldsymbol{\theta}^*\) the exact gradient vanishes, so the iteration stops whatever the sketch drew. Do you see why no such statement holds for sketch-and-solve? Pilanci & Wainwright prove that after \(T\) rounds the excess risk is small with high probability:
\[ L(\widehat{\boldsymbol{\theta}}_T) - L(\boldsymbol{\theta}^*) \le \chi^{2T}\,\|\mathbf{X}\boldsymbol{\theta}^*\|_2^2, \]
where the per-round contraction \(\chi\) is of order \(\sqrt{d/k}\). A sketch of size \(k \approx 6d\) contracts the error by a constant factor each round, so logarithmically many rounds finish the job. Every round releases a fresh sketch and an exact gradient. Both need a privacy account.
The Sketch Was Already a Mechanism
Now we look at what the sketch publishes. Each row of \(\mathbf{S}_t\mathbf{X}\) is \(\mathbf{X}^\top\mathbf{s}\) for an independent standard Gaussian \(\mathbf{s} \in \mathbb{R}^n\). A linear image of a standard Gaussian is Gaussian, with mean zero and covariance \(\mathbf{X}^\top\mathbf{I}_n\mathbf{X}\):
\[ \mathbf{X}^\top\mathbf{s} \sim \mathcal{N}(\mathbf{0}, \mathbf{X}^\top\mathbf{X}). \]
Releasing the sketch is releasing \(k\) i.i.d. samples from \(\mathcal{N}(\mathbf{0}, \mathbf{X}^\top\mathbf{X})\): the data enters only through its Gram matrix. Zeroing out person \(i\) dents that covariance by a rank-one outer product:
\[ \mathbf{X}^\top\mathbf{X} \;\longmapsto\; \mathbf{X}^\top\mathbf{X} - \mathbf{x}_i\mathbf{x}_i^\top. \]
How visible is a rank-one dent? In one dimension we can see the whole answer. There the sketch releases draws from \(\mathcal{N}(0, \lambda)\) with \(\lambda = \sum_i x_i^2\), and blanking person \(i\) shrinks the variance to \(\lambda - x_i^2\). An attacker has to tell those two Gaussians apart, and the only quantity in play is the person’s share \(x_i^2/\lambda\) of the variance: a share of one tenth moves the standard deviation by about five percent, and a share near zero leaves the two distributions almost indistinguishable. The general case measures the same share through the covariance.
Claim: fix a positive definite \(\boldsymbol{\Sigma} \in \mathbb{R}^{d \times d}\) and a vector \(\mathbf{x} \in \mathbb{R}^d\) whose share of the covariance, \(u = \mathbf{x}^\top\boldsymbol{\Sigma}^{-1}\mathbf{x}\), is below \(1\). The KL divergence of the dented Gaussian from the reference has a closed form:
\[ \mathrm{KL}\big(\mathcal{N}(\mathbf{0}, \boldsymbol{\Sigma} - \mathbf{x}\mathbf{x}^\top) \,\big\|\, \mathcal{N}(\mathbf{0}, \boldsymbol{\Sigma})\big) = \tfrac{1}{2}\big(-u - \ln(1 - u)\big), \]
which is \(\tfrac{u^2}{4} + O(u^3)\) for small \(u\). The linear term cancels, so a row’s visibility is second order in its share.
Proof
Start from the KL divergence between two centered Gaussians, for output covariance \(\boldsymbol{\Sigma}'\) against reference \(\boldsymbol{\Sigma}\):
\[ \begin{aligned} & \mathrm{KL}\big(\mathcal{N}(\mathbf{0}, \boldsymbol{\Sigma}') \,\big\|\, \mathcal{N}(\mathbf{0}, \boldsymbol{\Sigma})\big) \\ &= \tfrac{1}{2}\Big(\operatorname{tr}(\boldsymbol{\Sigma}^{-1}\boldsymbol{\Sigma}') - d + \ln\tfrac{\det \boldsymbol{\Sigma}}{\det \boldsymbol{\Sigma}'}\Big). \end{aligned} \]
Set \(\boldsymbol{\Sigma}' = \boldsymbol{\Sigma} - \mathbf{x}\mathbf{x}^\top\) and evaluate the two data-dependent pieces in turn. The trace uses the fact that the trace of an outer product is the matching inner product:
\[\begin{align} \operatorname{tr}(\boldsymbol{\Sigma}^{-1}\boldsymbol{\Sigma}') &= \operatorname{tr}(\mathbf{I}_d - \boldsymbol{\Sigma}^{-1}\mathbf{x}\mathbf{x}^\top) \\ &= d - \operatorname{tr}(\boldsymbol{\Sigma}^{-1}\mathbf{x}\mathbf{x}^\top) \\ &= d - \mathbf{x}^\top\boldsymbol{\Sigma}^{-1}\mathbf{x} \\ &= d - u. \end{align}\]
The determinant uses the matrix determinant lemma:
\[ \det(\boldsymbol{\Sigma} - \mathbf{x}\mathbf{x}^\top) = \det(\boldsymbol{\Sigma})\,(1 - \mathbf{x}^\top\boldsymbol{\Sigma}^{-1}\mathbf{x}) = \det(\boldsymbol{\Sigma})\,(1 - u), \]
so the log-determinant term is \(\ln\tfrac{\det \boldsymbol{\Sigma}}{\det \boldsymbol{\Sigma}'} = -\ln(1 - u)\). Substituting both pieces collapses the divergence:
\[ \tfrac{1}{2}\big((d - u) - d - \ln(1 - u)\big) = \tfrac{1}{2}\big(-u - \ln(1 - u)\big). \]
Expanding \(-\ln(1 - u) = u + \tfrac{u^2}{2} + \tfrac{u^3}{3} + \cdots\) cancels the linear term and leaves \(\tfrac{u^2}{4} + O(u^3)\).
The KL divergence is the mean of the privacy loss. Write \(\lambda_{\min}\) for the smallest eigenvalue of \(\mathbf{X}^\top\mathbf{X}\) and suppose it clears a floor \(\lambda_0 > 5/2\), the conditioning the analysis assumes. The Rayleigh quotient turns that floor into a bound on each row’s share, below the \(1\) the Claim requires:
\[ u = \mathbf{x}_i^\top(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{x}_i \le \frac{\|\mathbf{x}_i\|_2^2}{\lambda_{\min}} \le \frac{1}{\lambda_{\min}} \le \frac{1}{\lambda_0} < 1. \]
Each sketch row leaks a mean loss of about \(1/(4\lambda_0^2)\), the value \(\tfrac{u^2}{4}\) at \(u = 1/\lambda_0\), and the \(k\) independent rows add to \(k/(4\lambda_0^2)\). A mean is not yet an \((\varepsilon, \delta)\) guarantee. Sheffet [Sheffet, ICML ’17] discovered this route, and the tight Rényi analysis of Lev, Srinivasan, Shenfeld, Ligett, Sekhari & Wilson [LSSLSW, NeurIPS ’25] supplies the rest: the Rényi divergence of order \(\alpha\), whose \(\alpha \to 1\) limit is that mean, stays under the line \(\alpha k/(2\lambda_0^2)\) for every \(\alpha\) up to \(2\lambda_0/5\). The slope is twice the mean loss: the divergence bends upward, so a line pinned at the origin clears it only by starting out steeper. The standard conversion charges the slope, not the mean, plus twice the geometric mean of the slope and \(\ln(1/\delta)\):
\[ \varepsilon \le \frac{k}{2\lambda_0^2} + 2\sqrt{\frac{k}{2\lambda_0^2}\ln\tfrac{1}{\delta}} = \frac{k}{2\lambda_0^2} + \frac{\sqrt{2k\ln(1/\delta)}}{\lambda_0}. \]
The sketch alone is \((\varepsilon, \delta)\)-differentially private on a budget that shrinks as the floor grows: a person is a rank-one dent in a covariance, and a high floor hides that dent. Well-conditioned data is private on its own, with no noise added.
What if \(\lambda_{\min}\) falls short of the floor the budget demands? We add Gaussian noise to the sketch, \(\mathbf{S}\mathbf{X} + \eta\,\boldsymbol{\xi}\) with \(\boldsymbol{\xi} \in \mathbb{R}^{k \times d}\) i.i.d. standard Gaussian, and notice that this is the same mechanism run on a slightly larger dataset. We split a \(k \times (n+d)\) matrix \(\mathbf{S}'\) of i.i.d. standard Gaussians into its first \(n\) and last \(d\) columns:
\[ \mathbf{S}'\begin{bmatrix}\mathbf{X} \\ \eta \mathbf{I}_d\end{bmatrix} \;\overset{d}{=}\; \mathbf{S}\mathbf{X} + \eta\,\boldsymbol{\xi}. \]
Adding noise to the sketch is sketching the data with the \(d\) rows of \(\eta\mathbf{I}_d\) appended as phantom rows. They raise the Gram matrix to \(\mathbf{X}^\top\mathbf{X} + \eta^2\mathbf{I}_d\), so we pick the noise scale that lifts the smallest eigenvalue exactly to the floor:
\[ \eta^2 = \max\{\lambda_0 - \lambda_{\min},\, 0\}. \]
The noise fills the gap between the conditioning the data has and the conditioning privacy requires: it shrinks as \(\lambda_{\min}\) grows and vanishes at \(\lambda_{\min} = \lambda_0\). (\(\lambda_{\min}\) is data-dependent, so a small part of the budget goes to releasing it. Its sensitivity is at most \(1\) by Weyl’s inequality, and the release is deflated so the calibrated \(\eta\) is large enough with high probability.) That settles the curvature half. The gradient half is exact, and nothing exact is private on its own.
The Gradient Was Already a Sum
The exact gradient is a sum over people:
\[ \mathbf{X}^\top(\mathbf{y} - \mathbf{X}\widehat{\boldsymbol{\theta}}_t) = \sum_{i} \mathbf{x}_i \big(y_i - \mathbf{x}_i^\top \widehat{\boldsymbol{\theta}}_t\big). \]
A residual has no bound known in advance, so the sum has no sensitivity to calibrate. We clip each residual to \([-c, c]\). Then zeroing out person \(i\) moves the sum by a single term:
\[ \mathbf{x}_i \operatorname{clip}_c\big(y_i - \mathbf{x}_i^\top\widehat{\boldsymbol{\theta}}_t\big), \]
a vector of norm at most \(\|\mathbf{x}_i\|_2 \cdot c \le c\). A bounded-sensitivity sum is what the Gaussian mechanism handles: add \(\sigma \boldsymbol{\zeta}_t\) with \(\boldsymbol{\zeta}_t \sim \mathcal{N}(\mathbf{0}, \mathbf{I}_d)\) and \(\sigma\) calibrated to sensitivity \(c\) across the \(T\) rounds. This clip-and-noise treatment is how DP-SGD privatizes its gradients (Abadi, Chu, Goodfellow, McMahan, Mironov, Talwar & Zhang [AC+, CCS ’16]); here it lands on an algorithm whose curvature was already private for free.
One correction ties the two parts together. The phantom rows changed the problem: the augmented dataset \([\mathbf{X}; \eta\mathbf{I}_d]\) with responses \((\mathbf{y}; \mathbf{0})\) is ridge regression with penalty \(\eta^2\), whose gradient carries an extra \(-\eta^2\widehat{\boldsymbol{\theta}}_t\). Including that term keeps Newton’s gradient matched to Newton’s Hessian, so the fixed point sits at the ridge solution. (Dropping the term is also an option: the fixed point returns to \(\boldsymbol{\theta}^*\), though the contraction is slower, and the figure below does this. Privacy is untouched either way.)
The Algorithm
Both halves are mechanisms now, and one round assembles them from fresh randomness \(\mathbf{S}_t, \boldsymbol{\xi}_t, \boldsymbol{\zeta}_t\), a noisy sketch \(\widetilde{\mathbf{X}}_t \in \mathbb{R}^{k \times d}\) for the curvature and a noisy gradient \(\widetilde{\mathbf{G}}_t \in \mathbb{R}^d\) for the direction:
\[ \widetilde{\mathbf{X}}_t = \mathbf{S}_t \mathbf{X} + \eta\, \boldsymbol{\xi}_t, \qquad \widetilde{\mathbf{G}}_t = \mathbf{X}^\top \operatorname{clip}_c\big(\mathbf{y} - \mathbf{X}\widehat{\boldsymbol{\theta}}_t\big) - \eta^2 \widehat{\boldsymbol{\theta}}_t + \sigma \boldsymbol{\zeta}_t, \qquad \widehat{\boldsymbol{\theta}}_{t+1} = \widehat{\boldsymbol{\theta}}_t + \Big(\tfrac{1}{k}\widetilde{\mathbf{X}}_t^\top \widetilde{\mathbf{X}}_t\Big)^{-1} \widetilde{\mathbf{G}}_t. \]
In the diagram, the two private pathways meet in the update and loop back. Privacy is now pure accounting: the algorithm makes \(2T + 1\) private releases (one eigenvalue, \(T\) sketches, \(T\) gradients) whose Rényi curves add under composition. Iterating is nearly free: the noise scales grow like \(\sqrt{T}\) under composition while the optimization error shrinks like \(\chi^{2T}\), so the analysis needs only \(T\) of order \(\log(\max\{1, \lambda_{\min}/\lambda_0\})\). Extra rounds contract the error only while the data’s conditioning runs above the floor. Once any noise is being added, the analysis needs only one round. A constant \(T\) between 3 and 5 covers practice. One choice is still unexplained: the sketch touches \(\mathbf{X}\) and never \(\mathbf{y}\).
Why Sketching Only \(\mathbf{X}\) Matters
An earlier route, due to Sheffet and sharpened by [LSSLSW, NeurIPS ’25], sketches the concatenation \((\mathbf{X}, \mathbf{y})\) in one shot, so its noise is calibrated to the smallest eigenvalue of the joint Gram matrix \((\mathbf{X}, \mathbf{y})^\top(\mathbf{X}, \mathbf{y})\). That eigenvalue is dangerous. The smallest eigenvalue is at most the Rayleigh quotient at any unit vector, so we test the joint Gram matrix against the vector proportional to \((\boldsymbol{\theta}^*, -1)\):
\[\begin{align} & \lambda_{\min}\big((\mathbf{X}, \mathbf{y})^\top(\mathbf{X}, \mathbf{y})\big) \\ &\le \frac{\big\|(\mathbf{X}, \mathbf{y})\,(\boldsymbol{\theta}^*, -1)^\top\big\|_2^2}{\|(\boldsymbol{\theta}^*, -1)\|_2^2} \\ &= \frac{\|\mathbf{X}\boldsymbol{\theta}^* - \mathbf{y}\|_2^2}{1 + \|\boldsymbol{\theta}^*\|_2^2} \\ &= \frac{L(\boldsymbol{\theta}^*)}{1 + \|\boldsymbol{\theta}^*\|_2^2}. \end{align}\]
The better the linear fit, the smaller the residual \(L(\boldsymbol{\theta}^*)\), hence the smaller this joint eigenvalue and the more noise the joint sketch must add. A good fit is exactly the case the joint route handles worst. Sketching \(\mathbf{y}\) also puts the residual \(L(\boldsymbol{\theta}^*)\) into the error bound. No privacy budget removes that term.
The iterative Hessian sketch never touches \(\mathbf{y}\) with randomness. Pilanci & Wainwright kept the gradient exact because sketching it wrecks accuracy. A decade later, the same choice ties the noise to \(\lambda_{\min}\) of \(\mathbf{X}^\top\mathbf{X}\) alone and keeps the residual out of the error bound.
The Accounting
Here is the payoff [LSSLW, ICML ’26]. For the right \(k\), \(T\), and \(c\), IHM is \((\varepsilon, \delta)\)-differentially private, and with probability \(1 - \varrho\) its excess risk obeys a bound in a single multiplier:
\[ L(\widehat{\boldsymbol{\theta}}) - L(\boldsymbol{\theta}^*) = \widetilde{O}\Big(\gamma_{\mathrm{h}} \,(1 + \|\boldsymbol{\theta}^*\|_2^2)\, \min\{1, \gamma_{\mathrm{h}}/\lambda_{\min}\}\Big), \]
where \(\widetilde{O}\) hides logarithmic factors and the multiplier collects the privacy parameters:
\[ \gamma_{\mathrm{h}} = \frac{\sqrt{\max\{d,\, \ln(1/\varrho)\} \,\ln(1/\delta)}}{\varepsilon}. \]
Conditioning enters through the \(\min\): once \(\lambda_{\min}\) passes \(\gamma_{\mathrm{h}}\), the risk is quadratic in \(\gamma_{\mathrm{h}}\) rather than linear. AdaSSP satisfies a bound of the same shape with \(\gamma_{\mathrm{a}}\) in place of \(\gamma_{\mathrm{h}}\). (The max in \(\gamma_{\mathrm{h}}\) admits a third, data-dependent entry growing with \(\lambda_{\max} - \lambda_{\min}\); on typical data the other two dominate.) The whole improvement is the ratio of the two multipliers:
\[ \frac{\gamma_{\mathrm{a}}}{\gamma_{\mathrm{h}}} = \sqrt{\min\{d,\, \ln(1/\varrho)\}}. \]
Under AdaSSP the factors \(d\) and \(\ln(1/\varrho)\) multiply; under IHM only the larger of the two survives. That removes a multiplier as large as \(\sqrt{d}\) once \(\varrho\) is exponentially small.
The theory is visible in simulation. The setup is well conditioned: \(n = 2^{18}\) unit-norm rows spread over \(d = 128\) dimensions, so \(\lambda_{\min} \approx 1{,}960\), with \(\|\boldsymbol{\theta}^*\|_2^2 \approx 10\) and responses clipped to \([-1, 1]\). All three run at \(\delta = 1/n^2\), well under the \(1/n\) at which the slack could hide a whole released row, and IHM runs \(T = 5\) rounds. Every noise scale comes from the exact analytic Gaussian calibration of Balle & Wang [BW, ICML ’18] and the sketch accounting above, over 150 seeded trials.
In the plot, the one-shot sketch suffers both the joint-eigenvalue penalty and the residual amplification, and never recovers. AdaSSP and IHM both improve steadily, but IHM sits below AdaSSP across the entire range, by a factor of 12 at \(\varepsilon \approx 0.5\). Its calibrated sketch noise falls from \(\eta \approx 144\) at \(\varepsilon = 0.1\) to exactly zero past \(\varepsilon \approx 1.3\), where the curvature side privatizes itself and only the gradient noise remains. Across 33 real regression datasets [LSSLW, ICML ’26], IHM matches or beats AdaSSP and every sketching predecessor.
The lesson is not specific to regression. Randomized algorithms already inject calibrated randomness, and iterative ones already bound each person’s per-step influence. Those are the two raw materials of differential privacy. Sometimes the mechanism design is finished before the privacy question is asked, and all that remains is the accounting.