Iterative Hessian Mixing
13 min read
Suppose we want to publish a least squares fit trained on people’s data. An accurate fit is useful, but its coefficients can reveal information about the people whose records went into it. Differential privacy limits that disclosure by controlling how much the distribution of the released result can change when one person’s record changes. Adding enough randomness to protect each record can make the fit less accurate.
The iterative Hessian sketch already uses randomness to make least squares faster [PW, JMLR ’16]. Its random projection of the data can also provide privacy: when the data are well conditioned, changing one record has little effect on the projection’s distribution. The gradient used to update the fit still needs protection. We bound each person’s contribution to it by clipping, then add noise to the sum. Combining these two mechanisms gives iterative Hessian mixing [LSSLW, ICML ’26], a fast solver whose existing randomness helps meet the privacy requirement.
The Problem
What exactly has to stay hidden, and what does the standard method pay for hiding it? 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. IHM removes this factor. Removing it takes no new mechanism, only an account of what the sketch and the gradient each reveal about one person.
The Iterative Hessian Sketch
What is the iterative Hessian sketch, and why does a crude Hessian still converge? IHM starts from the following expansion of the least-squares 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 method 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 affects the convergence rate but not the fixed point: at \(\widehat{\boldsymbol{\theta}}_t = \boldsymbol{\theta}^*\) the exact gradient vanishes, so the iteration stops whatever the sketch drew. 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.
Privacy of the Sketch
How much does releasing the sketch \(\mathbf{S}_t\mathbf{X}\) reveal about one person? 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. \]
The Kullback–Leibler (KL) divergence quantifies the visibility of this rank-one change. In one dimension, the dependence is explicit. 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). \]
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 linear term cancels, so a row’s visibility is second order in its share: \(\tfrac{u^2}{4} + O(u^3)\) for small \(u\). The KL divergence is the mean of the privacy loss, while an \((\varepsilon, \delta)\) guarantee also requires a tail bound. 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)\). 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 uses the slope, rather than 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.
If \(\lambda_{\min}\) falls short of the required floor, 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. 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.)
Privatizing the Gradient
What about the exact gradient, which the sketch leaves untouched? 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.
The phantom rows change 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
What does one round look like with both parts private, and what does iterating cost? One round combines 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. 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.
Sketching Only \(\mathbf{X}\)
Why sketch only \(\mathbf{X}\) and not the responses along with it? 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})\). The smallest eigenvalue is at most the Rayleigh quotient at any unit vector, so 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}\]
A smaller residual \(L(\boldsymbol{\theta}^*)\) makes this joint eigenvalue smaller and forces the joint sketch to add more noise. Sketching \(\mathbf{y}\) also puts the residual into the error bound, where no privacy budget removes it.
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
How much better is the guarantee than AdaSSP’s? Lev, Shenfeld, Srinivasan, Ligett & Wilson prove that for the right \(k\), \(T\), and \(c\), IHM is \((\varepsilon, \delta)\)-differentially private [LSSLW, ICML ’26]. 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 the gap between the largest and smallest eigenvalues; on typical data the other two dominate.) The ratio of the two multipliers is:
\[ \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.
IHM privatizes an existing randomized iterative algorithm by calibrating its randomness to each person’s per-step influence and composing the resulting privacy losses.