Iterative Hessian Mixing
13 min read
We want to fit a model to people’s data without letting the fitted model reveal too much about any one person. Differential privacy makes this requirement precise: the distribution of the released result must change little when one person’s record is removed. Adding enough noise to hide that influence can make the fit less accurate.
Some solvers already use randomness to make fitting faster. Surprisingly, iterative Hessian mixing (IHM) can use that same randomness to help satisfy the privacy requirement [LSSLW, ICML ’26].
The Problem
Consider fitting a linear model by least squares, which minimizes the sum of squared prediction errors. The dataset has one row per person: features \(\mathbf{x}_i \in \mathbb{R}^d\) and a response \(y_i\). Collect them in a matrix \(\mathbf{X} \in \mathbb{R}^{n \times d}\) and a vector \(\mathbf{y} \in \mathbb{R}^n\), where \(n\) is the number of people. 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}^*)\).
Normalize each row so that \(\|\mathbf{x}_i\|_2 \le 1\) and \(|y_i| \le 1\). 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. IHM handles the curvature and the gradient separately, using the random sketch to protect the curvature. We can compare its error with AdaSSP once those two privacy costs are clear.
The Iterative Hessian Sketch
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 and keeps the gradient exact [PW, JMLR ’16]. It 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. They need different privacy treatments:
| Quantity | What the fast solver computes | What privacy requires |
|---|---|---|
| Curvature | A random sketch of the Gram matrix | Account for the sketch’s randomness; add noise if needed |
| Gradient | The exact gradient | Clip each person’s contribution and add noise |
The sketch’s existing randomness can help protect the curvature, while the exact gradient still needs a separate adjustment.
Privacy of the Sketch
To account for privacy, we need the distribution of the released sketch. 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
The sketch leaves the gradient exact, so that part of the computation still needs protection. 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
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}\)
Keeping the responses out of the sketch separates the cost of privacy from how well the model fits. 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
We can now compare the error guarantee with the AdaSSP baseline. AdaSSP’s 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.
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.