Revisiting inverse Hessian vector products for calculating influence functions

## Abstract

Abstract Influence functions are a popular tool for attributing a model’s output to training data. The traditional approach relies on the calculation of inverse Hessian-vector products (iHVP), but the classical solver “Linear time Stochastic Second-order Algorithm” (LiSSA, Agarwal et al. [ 2017 ] ) is often deemed impractical for large models due to expensive computation and hyperparameter tuning. We show that the three hyperparameters — the scaling factor, the batch size, and the number of steps — can be chosen depending on the spectral properties of the Hessian, particularly its trace and largest eigenvalue. By evaluating with random sketching [Swartworth and Woodruff, 2023 ] , we find that the batch size has to be sufficiently large for LiSSA to converge; however, for all of the models we consider, the requirement is mild. We confirm our findings empirically by comparing to Proximal Bregman Retraining Functions (PBRF, Bae et al. [ 2022 ] ). Finally, we discuss what role the inverse Hessian plays in calculating the influence.

## 1 Introduction

Deep neural networks have seen many impressive results in the past years, but researchers and practitioners have little understanding what happens inside the models and how they learn to predict. While there are exist various methods to interpret the internal computations in an understandable to a human way, influence functions attempt to explain model behaviour by attributing model predictions (or generations) to particular examples in the training data.

Koh and Liang [] introduce Hessian-based influence functions in order to approximate the effect of removal of one training point from the training set, which we refer to as leave-one-out retraining. The formula for influence calculation is derived from the second-order Taylor approximation of the loss, hence the Hessian and the gradient of the training point are sufficient for calculation. Koh and Liang [] demonstrate various applications of influence functions such as explaining of model outputs through data attribution, repairing mislabeled data, and backdoor attacks.

Basu et al. [] criticize influence functions for poor approximation of leave-one-out retraining as depth and width of neural networks increase. As a solution, Bae et al. [] propose two fixes: replacing the Hessian (that possibly has negative eigenvalues) with well-behaved Gauss-Newton Hessian [Martens, ], and replacing the leave-one-out retraining (which itself is not a well-defined objective) with Proximal Bregman Retraining Functions (PBRF). They demonstrate that the latter do not suffer from the randomness introduced by model initialization and data sampling, and they argue can serve as gold standard when evaluating influence function approximation methods. In this paper, we focus on this particular formulation of influence functions, where the Hessian is replaced with Gauss-Newton Hessian, and the PBRF serves as the ground truth in validation.

The calculation of influence functions, as introduced by Koh and Liang [], requires approximation of inverse Hessian-vector products. Given the dimension of modern deep models and the size of the training dataset, it can be a hard problem. As an alternative to traditional conjugate gradient method, Koh and Liang [] proposed a stochastic iterative approach called “Linear time Stochastic Second-Order Algorithm”, or in short LiSSA [Agarwal et al., ]. This algorithm requires calculating a sampled Hessian-vector product at each iteration, where in the formulation of Koh and Liang [], the batch size per sample can be as little as just one training point. There are two additional hyperparameters involved — the scaling factor and the number of steps in LiSSA, with little direction of how to choose them in practice. Basu et al. [] criticize LiSSA, suggesting that it lacks convergence for deep networks with large number of parameters. In particular, the method deemed impractical due to the need for expensive hyperparameter search. In this paper, we carefully analyze the convergence of LiSSA and find that the choice of all three hyperparameters, including the batch size, depends on the properties of the Gauss-Newton Hessian, namely, its trace and largest eigenvalue. Since the size of the Hessian is very large (number of parameters to square), we evaluate these two statistics with *random sketching* [Swartworth and Woodruff, ], which only requires estimation of Hessian-vector products in the process. We report these statistics and the corresponding requirements for some open sourced vision and language models. We also note that, contrary to a common belief, we find that the batch size has to be sufficiently large for the algorithm to converge. However, this requirement is mild, and for language models it is particularly redundant.

Some attempts to avoid calculating inverse Hessian vector products have been made in the literature. Schioppa et al. [] suggest to truncate the spectrum of the Hessian, see also Fisher et al. [], Grosse et al. [].
When it comes to language models, most of the recent literature is using gradient-based influence functions [Xia et al., , He et al., , Chhabra et al., ]. These are typically focused on the finetunning stage, and often this choice is motivated by simpler and faster implementation.
In Section 5, we highlight the difference between the results of these two ways of calculating influence functions.

Almost exclusively, Grosse et al. [] calculate Hessian-based influence of pretraining data for LLMs. Their analysis is restricted to MLPs of the transformer and they impose a block-wise structure onto the Hessian. Although we do not advocate against such structural assumptions, our work suggests that running the plain and model-agnostic LiSSA can be feasible, given that we avoid the hyperparameter search.
In our implementation, we follow the Hessian-free approach of Martens et al. [] using finite differences, where an in-batch GNH-vector product is calculated with three forward propagations and one backward propagation. For example, our code allows to run LiSSA for 7B language models on a 4$\times$A100 GPU node.
We share the code here [https://github.com/yklochkov-bytedance/gnhtools](https://github.com/yklochkov-bytedance/gnhtools).

## 2 Background and notation

Influence functions [Koh and Liang, ] are calculated under the assumption that the optimized parameter $\theta$ of the model delivers minimum to the training loss,

$$ $\theta^{*}=\arg\min_{\theta}\frac{1}{|\mathcal{D}_{tr}|}\sum_{(\mathbf{x},y) \in\mathcal{D}_{tr}}\ell(\mathbf{x},y;\theta),$ (1) $$

where for classification tasks, $(\mathbf{x},y)$ are input and label pair, and for language modeling tasks, consists of context and next word token. That is, given a sequence $s=(s_{1},\dots,s_{l})$, the dataset $\mathcal{D}_{tr}$ consists of pairs $\mathbf{x}=(s_{1},\dots,s_{t-1})$ and $y=s_{t}$. Let us fix a point $(\mathbf{x}_{m},y_{m})\in\mathcal{D}_{tr}$, and for a small perturbation weight $\epsilon>0$ consider

$$ $\theta^{*}(\epsilon)=\arg\min_{\theta}\frac{1}{|\mathcal{D}_{tr}|}\sum_{( \mathbf{x},y)\in\mathcal{D}_{tr}}\ell(\mathbf{x},y;\theta)+\epsilon\,\ell( \mathbf{x}_{m},y_{m};\theta).$ (2) $$

Then, the *influence of a training point* $(\mathbf{x}_{m},y_{m})$ on the parameter is denoted as

$$ $\displaystyle\mathcal{I}((\mathbf{x}_{m},y_{m}))$ $\displaystyle=\frac{d\theta^{*}(\epsilon)}{d\epsilon}\Big{|}_{\epsilon=0}$ $\displaystyle=-H^{-1}\nabla\ell(\mathbf{x}_{m},y_{m};\theta^{*}),$ (3) $$

where $H$ denotes the *population Hessian*, that is

$$ $H=\frac{1}{|\mathcal{D}_{tr}|}\sum_{(\mathbf{x},y)\in\mathcal{D}_{tr}}\nabla^{ 2}\ell(\mathbf{x},y;\theta^{*})\,.$ $$

Furthermore, suppose we have a set of predictions $z_{test}=(\mathbf{x}_{test},\hat{y}_{test})$, and let $f((\mathbf{x},y),\theta)=\log p(y|\mathbf{x};\theta)$ are the log probability according to the trained model. Then, the influence of a *training point $z_{train}=(\mathbf{x}_{m},y_{m})$ on a prediction $z_{test}$* is denoted as

$$ $\mathcal{I}(z_{train},z_{test})=-\nabla f(z_{test},\theta)^{\top}H^{-1}\nabla \ell(\mathbf{x}_{m},y_{m};\theta^{*})$ $$

For language models, we calculate the influence for completion following Grosse et al. []. Let $s=(s_{1},\dots,s_{p})$ be a prompt and $\hat{s}=(\hat{s}_{1},\dots,\hat{s}_{c})$ be a completion. Then, we calculate the influence for average log-probability of predicted tokens

$$ $f(s,\hat{s};\theta)=\frac{1}{c}\sum_{j=1}^{c}\log p(\hat{s}_{j}|s_{1}\dots s_{ p}\hat{s}_{1}\dots\hat{s}_{j-1};\theta)\,.$ $$

#### PBRF as ground truth.

The inverse problem $H^{-1}\nabla\ell(z_{train};\theta)$ can be difficult to perform due to degenerate eigenvalues of $H$. Koh and Liang [] propose to use a damping parameter $\lambda>0$ and instead invert a regularized matrix $(H+\lambda I)^{-1}$.
However, such matrix can still be degenerate due to possibly negative eigenvalues of the Hessian of a non-convex loss, which are indeed observed in practice and are not necessarily small [Sagun et al., , Schioppa, ].
Motivated by natural descent methods, Bae et al. [] propose to replace it with *Gauss-Newton Hessian* (GNH), which is denoted as follows. Suppose that the loss has the form

$$ $\ell((\mathbf{x},y);\theta)=\ell(h(\mathbf{x};\theta),y),\qquad\ell(h,y)=-\log (\mathrm{sf}(h)_{y}),$ $$

where $h(\mathbf{x};\theta)\in\mathbb{R}^{K}$ is the logit function and $\mathrm{sf}(h)_{j}=\exp(h_{j})/\left(\sum_{j=1}^{K}\exp(h_{j})\right)$ is the standard softmax function. Then, the GNH has the form

$$ $H=\frac{1}{|\mathcal{D}_{tr}|}\sum_{(\mathbf{x},y)\in\mathcal{D}_{tr}}[J_{ \theta}h(\mathbf{x};\theta)]\nabla_{h}^{2}\ell(h(\mathbf{x};\theta),y)[J_{ \theta}h(\mathbf{x};\theta)]^{\top}$ (4) $$

For the Cross-Entropy loss, we have the identity $\nabla_{h}^{2}\ell(h(\mathbf{x};\theta),y)=\mathrm{Diag}(\mathrm{sf}(h))-
\mathrm{sf}(h)\mathrm{sf}(h)^{\top}$.

Furthermore, Bae et al. [] show that if the Hessian is replaced with the Gauss-Newton Hessian, the influence functions (3) approximate a different retraining functions called Proximal Bregman Retraining Functions (PBRF). These correspond to retraining of the Proximal Bregman Objective (PBO) on training point $(\mathbf{x}_{m},y_{m})$ reads as follows,

$$ $\theta_{PBRF}(\epsilon)=\arg\min_{\theta}\frac{1}{|\mathcal{D}_{tr}|}\sum_{( \mathbf{x},y)\in\mathcal{D}_{tr}}D(h(\mathbf{x};\theta),h(\mathbf{x};\theta^{ \star}),y)+\epsilon\ell((\mathbf{x}_{m},y_{m});\theta)+\frac{\lambda}{2}\| \theta-\theta^{\star}\|^{2},$ (5) $$

where $D(h,h^{\prime},y)=\ell(h,y)-\ell(h^{\prime},y)-(h-h^{\prime})^{\top}\nabla_{h}
\ell(h^{\prime},y)$ is the Bregman divergence. Comparing PBO with the objective in (2), we see that the proximity penalty $\frac{\lambda}{2}\|\theta-\theta^{\star}\|^{2}$ takes into account the damping parameter, while replacing the loss with the Bregman divergence accounts for potential lack of convergence, i.e. we no longer need to assume that the training of the original parameter converges to global minimum of the loss as in (1). In addition, Bae et al. [] find that PBRFs are a more reliable objective compared to traditional retraining, which is known to produce different outputs [Basu et al., ].

Bae et al. [] argue that PBRF is suitable ground truth objective for validation of influence function estimation algorithms. For instance, they use it for empirical confirmation of their ad-hoc algorithm in Grosse et al. []. Following them, we refer to PBRF as a ground truth influence in order to confirm our findings empirically.

#### Iterative inverse Hessian-vector products.

For calculating these inverse Hessian-vector products (iHVP) of form $\mathbf{u}=(H+\lambda)^{-1}\mathbf{g}$, Koh and Liang [] propose to use a variant of Linear time Stochastic Second-Order Algorithm (LiSSA, Agarwal et al. []), that consist of the iterations

$$ $\mathbf{u}^{t}=\mathbf{g}+(I-\eta(\tilde{H}^{t}+\lambda))\mathbf{u}^{t-1}, \qquad t=1,\dots,T$ (6) $$

where $\tilde{H}^{t}$ is an in-batch estimate of $H$. Ideally, the scaling parameter $\eta>0$ needs to be chosen to ensure that $\eta H$ is a contraction, however, it requires knowning the larges eigenvalue of $H$. To this day, LiSSA is often discarded due to it’s hyperparameters, which are not trivial to tune when one does not have a clear objective Basu et al. []. In particular, the three hyperparameters — $\eta$, $T$, and the batch size, are often chosen without any directive, and the resulting estimate is deemed unreliable Basu et al. [].

The LiSSA updates (6) are equivalent to stochastic gradient descent (SGD) with step size $\eta$ for the quadratic objective $\min_{\mathbf{u}}\frac{1}{2}\mathbf{u}^{\top}(H+\lambda)\mathbf{u}-\mathbf{u}^
{\top}\mathbf{g}$ [Fisher et al., ]. In theory, mini-batch SGD is known to work at least as well as full gradient, even in terms of number of updates [Harvey et al., ]. However, larger batch sizes are often preferred by practitioners. The problem stems from theoretical results relying on a notion of uniform smoothness, which is not observed in practice, where an in-batch function is usually not as smooth as the average over the whole dataset [Tang et al., ]. Furthermore, optimal choice of step size $\eta$ and the number of steps $T$ depend on the largest eigenvalue of the Hessian $\lambda_{\max}(H)$. Rather than conducting hyperparameter tuning, we suggest to evaluate the largest eigenvalue $\lambda_{\max}(H)$ directly, which allows us to run the LiSSA just once.

#### Hessian-vector products.

The updates (6) involve calculation of in-batch Hessian-vector products $\tilde{H}_{t}\mathbf{u}$. Expanding the expression for GNH (4), we observe that for a batch of data $B=\{(\mathbf{x},y)\}$, we have

$$ $\tilde{H}_{t}\mathbf{u}=\frac{1}{|B|}\sum_{(\mathbf{x},y)\in B}[J_{\theta}h( \mathbf{x};\theta)]\nabla_{h}^{2}\ell(h(\mathbf{x};\theta),y)[J_{\theta}h( \mathbf{x};\theta)]^{\top}\mathbf{u}$ (7) $$

Here $[J_{\theta}h(\mathbf{x};\theta)]^{\top}\mathbf{u}$ is a directional derivative of vector-function $h(\mathbf{x};\theta)$. Calculating the directional derivatives per each example in the batch precisely may be prohibitively expensive. Instead, we suggest to approximate it by finite differences,

$$ $[J_{\theta}h(\mathbf{x};\theta)]^{\top}\mathbf{u}=\frac{d}{d\delta}h(\mathbf{x };\theta+\delta\mathbf{u})\Big{|}_{\delta=0}\approx\frac{h(\mathbf{x};\theta+ \delta\mathbf{u})-h(\mathbf{x};\theta-\delta\mathbf{u})}{2\delta},$ $$

where $\delta$ is a small value, which we fix to $\delta=0.01$ in our experiments. Then, we can approximate the in-batch GNH-vector prodcut by using three forward propagations and one backward propagation:

$$ $\tilde{H}_{t}\mathbf{u}\approx\nabla_{\theta}\frac{1}{|B|}\sum_{(\mathbf{x},y) \in B}h(\mathbf{x};\theta)^{\top}S_{h}\left\{\frac{h(\mathbf{x};\dot{\theta}+ \delta\mathbf{u})-h(\mathbf{x};\dot{\theta}-\delta\mathbf{u})}{2\delta}\right\},$ $$

where $\dot{\theta}$ indicates that we do not calculate the derivative through this parameter, and $S_{h}=\mathrm{Diag}(\mathrm{sf}({h}))-\mathrm{sf}({h})\mathrm{sf}({h})^{\top}$ is also fixed. Thus, we simply backpropogate through a weighted sum of logits in the batch $h(\mathbf{x};\theta)$, with weights depending on the matrices $S_{h}$ and finite differences $(h(\mathbf{x};\dot{\theta}+\delta\mathbf{u})-h(\mathbf{x};\dot{\theta}-\delta
\mathbf{u}))/(2\delta)$. The latter two can be calculated in a gradient free manner with three forward propagations.
We note that incorporating finite differences for Hessian-vector products calculation has previously been done in Martens et al. [], Martens and Sutskever [].

## 3 Convergence of LiSSA and choice of hyperparameters

In order to carefully analyze the convergence of LiSSA iterations (6) we reformulate them as SGD updates. Observe that the result of iHVP applied to a gradient $\mathbf{g}$, $\mathbf{u}^{\star}=(H+\lambda)^{-1}\mathbf{g}$, delivers minimum to the following objective

$$ $\mathbf{u}^{\star}=\arg\min_{\mathbf{u}}L(\mathbf{u}),\qquad L(\mathbf{u}):= \frac{1}{2}\mathbf{u}^{\top}{H}\mathbf{u}+\frac{\lambda}{2}\|\mathbf{u}\|^{2}- \mathbf{u}^{\top}\mathbf{g}\,.$ (8) $$

With appropriate scaling, the LiSSA updates are equivalent to SGD with step size $\eta$ and the gradient calculated on in-batch loss $\tilde{L}_{t}(\mathbf{u})$, where the Gauss-Newton Hessian $H$ is replaced with unbiased estimate $\tilde{H}_{t}$ calculated over a random batch, turning the updates in (6) into

$$ $\mathbf{u}^{t}=\mathbf{u}^{t-1}-\eta\left[(\tilde{H}^{t}+\lambda)\mathbf{u}^{t -1}-\mathbf{g}\right],\qquad t=1,\dots,T$ $$

where the scaling parameter in (6) now plays the role of a learning rate.
Convergence of SGD is well studied in the literature, with the recommending step size $\eta$ typically depending on the smoothness of $L(\mathbf{u})$, which in our case equals to $\lambda_{\max}(H)$ Bubeck et al. []. Although in theory, mini-batch SGD is generally considered more efficient than full-batch, in practice mini-batch SGD often performs poorly due to the difference in smoothness of the population objective $L(\mathbf{u})$ and the in-batch objective $\tilde{L}(\mathbf{u})$ as pointed out by Tang et al. []. Instead, they formulate their bounds in terms of expected smoothness of $\tilde{L}_{t}(\mathbf{u})$. We derive the following results that takes into account this difference for the quadratic optimization (8).

###### Theorem 1 .

Suppose, $\eta<1/(\lambda_{\max}(H)+\lambda)$. Then, we have convergence in-expectation

$$ $\|\mathbf{E}\mathbf{u}^{t}-\mathbf{u}^{\star}\|\leq(1-\lambda\eta)^{t}\| \mathbf{u}^{0}-\mathbf{u}^{\star}\|\,.$ $$

Furthermore, assume that $\eta>0$, $\delta\in(0,1)$ are such that

$$ $(1-\eta(H+\lambda))^{2}+\eta^{2}(\mathbf{E}\tilde{H}_{t}^{2}-H^{2})\preceq(1- \delta)I.$ (9) $$

Then,

$$ $\mathbf{E}\|\mathbf{u}^{t}-\mathbf{u}^{\star}\|^{2}\leq(1-\delta)^{t}\left(2\| \mathbf{u}^{0}-\mathbf{u}^{\star}\|^{2}+\|\mathbf{u}^{\star}\|^{2}\right)+ \delta^{-1}\eta^{2}\tilde{\Delta},$ $$

where we interpret $\tilde{\Delta}=\mathbf{E}\|(H-\tilde{H}_{t})\mathbf{u}^{*}\|^{2}$ as a sampling error.

We conclude that although the updates of LiSSA are designed to be unbiased, it is not always guaranteed to converge even if $\eta<1/\lambda_{\max}(H)$.
In Section B.4, we show a counter-example where the difference $\mathbf{E}\|\mathbf{u}^{t}-\mathbf{u}^{\star}\|^{2}$ is not guaranteed to converge whenever the requirement (9) does not hold.
The principal difference comes from the matrix $\mathbf{E}\tilde{H}_{t}^{2}-H^{2}$, which can be interpreted the sampling gap, i.e. the larger the batch size, the closer the sampled squared Hessian is to the population squared Hessian. In particular, as we increase the batch size $|B|\rightarrow\infty$, sooner or later we expect to have that $\mathbf{E}\tilde{H}_{t}^{2}\approx H^{2}$. Thus, the batch size has direct impact on the convergence of LiSSA. Although it is hard to assess the inequality (9) directly, we find that in many cases we can rely on the following simple condition

$$ $\boxed{\mathbf{E}\tilde{H}_{t}^{2}-H^{2}\preceq\frac{C}{|B|}\mathrm{Tr}(H)H,}$ (C.1) $$

where $|B|$ is the batch size, which for language models corresponds to the total in-batch *number of tokens*. For simplicity, we assume this number to be the same in each batch $B$. We think of $C$ as a constant moderately larger than $1$, e.g. $C=2$.

In particular, we find this condition to hold for classification with independent sampling (assuming certain properties of the gradients’ distribution). For the case of language modeling, the batch consists of tokens sampled per sequence. We show that the condition holds if the gradients corresponding to different tokens within the same sequence have (mostly) little correlation. This is reasonable to expect due to the large dimension of the parameter, and we also note that a similar assumption appears in Tang et al. [] in the context of imaging inverse problems.
In addition, we provide a simple empirical test where we compare the traces of LHS and RHS in (C.1), and it confirms the inverse scaling with batch size for both classification and language modeling tasks. See Section B.3 in the appendix for these derivations and experiments.

Under condition C.1, we can rewrite Theorem 1 in a simplified form with an exact requirement for a sufficiently large batch size.

###### Corollary 1 .

Suppose that C.1 holds. Let us choose the hyperparameters

$$ $\eta=1/(\lambda_{\max}(H)+\lambda),\qquad|B|\geq C\mathrm{Tr}(H)/\lambda_{\max }(H)$ $$

Then,

$$ $\mathbf{E}\|\mathbf{u}^{t}-\mathbf{u}^{*}\|^{2}\lesssim(1-\lambda\eta)^{2t}(\| \mathbf{u}^{0}\|^{2}+\|\mathbf{u}^{\star}\|^{2})+\frac{\eta^{2}\mathrm{Tr}(H)} {|B|}\mathbf{g}^{\top}(H+\lambda)^{-1}\mathbf{g}$ (10) $$

and the algorithm converges in $T=\Omega(1/(\eta\lambda))$ steps.

In our error bound above, the first term depends on the learning rate and the number of steps, and we can say that it measures how quickly we converged to the solution, therefore we label it *convergence error*. The second term depends directly on the batch size and we label it *sampling error*. It does not depend on the number of steps performed and comes from the difference between the sampled HVP and population HVP. In particular, it is trivial to see that it corresponds to the variance of a single update, $\mathbf{E}\|\mathbf{u}^{t}-\mathbf{E}[\mathbf{u}^{t}|\mathbf{u}^{t-1}]\|^{2}$, in the limit $\mathbf{u}^{t}\rightarrow\mathbf{u}^{\star}$. Notice that although the *convergence error* does not depend on the batch size explicitly, we have to satisfy the condition $|B|\geq C\mathrm{Tr}(H)/\lambda_{\max}(H)$ in order to converge. We also note that in terms of batch size requirement our result conforms with Tang et al. []. However, they only make a qualitative characterization that mini-batch SGD is applicable in problems where $H$ has fast decaying eigenspectrum.

###### Remark 1 .

Notice that Theorem 1 can hold for arbitrary matrix $H$ and its sampling counterpart, such as the original Hessian $\mathbf{E}\nabla^{2}\ell(z;\theta)$, and arbitrary unbiased estimates $\tilde{H}_{t}$, assuming the batches are drawn independently of each other and from the same distribution. However, condition (C.1) and, correspondingly, Corollary 1 are only expected to hold for the Gauss-Newton Hessian.

**Table 1: Gauss-Newton-Hessian statistics and recommended hyperparameters. Statistics are calculated on ImageNet (IN) on vision and Open-Web-Text-2 (OWT) on language models. The arrow $\uparrow$ indicates a lower bound, and $\downarrow$ indicates that it is an upper bound.**
|  |  |  |  | $\lambda_{\max}(H)$ | recommended<br>hyperparameters^* |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- |
| Model | Size | Data | $\tfrac{1}{N}\mathrm{Tr}(H)$ | from<br>sketching | $\downarrow\eta$ | $\uparrow|B|$ | $\uparrow T$ |
| ResNet-18 | 11M | IN | $(1.32\pm 0.00)\times 10^{-3}$ | $\approx 270$ | 0.003 | 100 | 150 |
| ResNet-50 | 25M | IN | $(8.17\pm 0.11)\times 10^{-4}$ | $\approx 470$ | 0.002 | 5 | 200 |
| OPT | 1.3B | OWT | $(9.28\pm 0.35)\times 10^{-6}$ | $\approx 780$ | 0.001 | 30 | 500 |
| Llama-1 | 7B | OWT | $(5.69\pm 0.67)\times 10^{-6}$ | $\approx 1600$ | 0.0005 | 50 | 1000 |
| Mistral | 7B | OWT | $(8.18\pm 0.13)\times 10^{-5}$ | $\approx 5600$ | 0.0002 | 200 | 2000 |

###### Theorem 1 .

###### Corollary 1 .

###### Remark 1 .

#### Empirical analysis of eigenvalue statistics.

It is apparent that to choose hyperparameters correctly, we need to evaluate the statistics $\lambda_{\max}(H)$ and $\mathrm{Tr}(H)$. Since there is no way to calculate the Hessian explicitly, we resort to random feature methods that only require evaluation of HVPs.

Evaluating the trace is straightforward. We generate a series of quadratic forms $(\mathbf{g}_{i}^{\top}\tilde{H}_{i}\mathbf{g}_{i})_{i=1}^{N}$, which consists of calculating HVPs and dot products for Gaussian vectors $\mathbf{g}\sim\mathcal{N}(0,\tfrac{1}{N}I)$. Their mean estimates the trace $\mathbf{E}\mathbf{g}^{\top}\tilde{H}\mathbf{g}=\mathbf{E}\mathbf{g}^{\top}H
\mathbf{g}=\frac{1}{N}\mathrm{Tr}(H)$. Due to independence of observations, we can also evaluate standard error of this estimator.
For evaluating the largest eigenvalue we use *random sketching*. That is, we evaluate the matrix $\hat{H}=\Phi H\Phi^{\top}$, with $\Phi\in\mathbb{R}^{d\times N}$ generated in a way such that $\Phi_{i,j}\sim\mathcal{N}(0,1/d)$. It is known, such sketches can evaluate the top eigenvalues of the original matrix $\hat{\lambda}_{i}(H)=\lambda_{i}(\hat{H})-\frac{1}{d}\mathrm{Tr}(\hat{H})$ [Swartworth and Woodruff, ], with error of estimation negligible for the top eigenvalue. See calculation details in the appendix, Section C.

We report the results of these evaluations in Table 1 for 2 ResNets and 3 open-sourced language models. We also show recommendations for the choice of hyperparameters based on Corollary 1. Notice that contrary to the original idea of SGD, in all cases the recommended batch size is larger than $1$.
However, recall that for language modeling, the batch size is the amount of tokens in a batch, and the recommended values are smaller than a typical context length. Thus, the LiSSA can work with just one sequence per batch. We also note that the recommendation is only a lower bound that ensures that LiSSA does not diverge, and increasing the batch size further makes the sampling error smaller (second term in (10)).

## 4 Empirical validation

Figure: Figure 1: Comparison of PBRF and LiSSA influence. The first row shows examples of training images. Below, the $x$-axis represents LiSSA influences, and the $y$-axis represents the PBRF influences corresponding to each training image and 500 test images. The second row is for ResNet-18, and the third row is for ResNet-50.
Refer to caption: extracted/5879605/figures/pbrf_comparison.png

Figure: Figure 2: Convergence of LiSSA for ResNet-18 with different batch size configurations. We calculate the correlation between test influences at steps 1..1000 of LiSSA. The result for the small batch size of 10 is averaged over 10 trials, so that the amount of data used in the middle and rightmost figures is the same.
Refer to caption: extracted/5879605/figures/conv.png

We now conduct empirical validation of our theoretical results. In particular, we want to check two things. Firstly, we want to demonstrate that when the parameters are chosen according to Table 1, the LiSSA converges as expected. For ground truth, we calculate PBRF for selected training examples Bae et al. []. Secondly, we want to empirically confirm that the requirement on sufficiently large batch size is indeed important.

We compare LiSSA and PBRF for the ResNet-18 and ResNet-50 models and randomly selected 25 training and 500 test images from the ImageNet dataset. For each training image we calculate iHVP $s_{train}=(H+\lambda)^{-1}\nabla\ell(x_{train},y_{train})$ using LiSSA with hyperparameters from Table 1, and calculate the 500 influences $s_{train}^{\top}\nabla\ell(x_{test},y_{test})$. We also calculate the PBRF by finetuning the model with SGD on Proximal Bregman objective. For PBRF, we match the batch size, number of steps, and the learning rate to LiSSA. We take $\epsilon=$1e-8 in (5) and optimize the PBO using double precision to avoid float overflow. We show results for ResNet-18 and ResNet-50 for selected 5 images, and the full list is shown in the appendix, Section E. We observe three cases: 1) LiSSA approximates PBRF, i.e. scatter plot concentrates along the dashed line $x=y$; 2) both LiSSA and PBRF have very low values poorly distinguishable from zero; 3) both LiSSA and PBRF have high value and do not approximate each other. In the latter case, we can argue that the PBO finetunning stirs away the model too far for the quadratic approximation to hold.

Furthermore, we confirm that the batch size matters not only for the sampling error in (10), but also for the speed of convergence. Let us take ResNet-18 with damping parameter and run the LiSSA algorithm for 1000 steps. According to Table 1, the recommended batch size is equal to 100. We suggest to consider three set-ups: batch size 2x larger than recommended, batch size is equal to 100, and batch size 10x smaller than recommended, with the result averaged over 10 independent runs. In the latter and the former cases, the result is obtained through iterating in the same amount of data, which allows to equalize the sampling error of a single update. Figure 2 reports the correlation between influences calculated for 500 test images at different steps. As we can see the correlation converges to 1 faster for the two cases where the batch size is greater or equal to recommended 100, despite averaging over 10 trials for the small batch.

## 5 What is the role of inverted Hessian?

In the context of language models, the focus in the current literature is mostly on gradient-based influence functions [Xia et al., , He et al., , Chhabra et al., ]. Often this choice is motivated by simpler and faster implementation. Due to the high cost of Hessian-based influence calculation, it is natural to ask what are the benefits compared to the gradient-based influence. We conduct a simple experiment in an attempt to understand what is left out of the consideration when relying only on gradient dot products.

Consider the eigenvalue decomposition of the Gauss-Newton Hessian $H$,

$$ $H=\sum_{j=1}^{N}\lambda_{j}v_{j}v_{j}^{\top},$ $$

where $v_{j}$ are orthogonal and normalized and $Hv_{j}=\lambda_{j}v_{j}$. If we represent a gradient $g=\nabla\ell(z_{test})$ in this eigenbasis, the iHVP simply reweights the coefficient according to how large the eigenvalue is,

$$ $g=\sum_{j}\langle g,v_{j}\rangle v_{j},\qquad\lambda(H+\lambda)^{-1}g=\sum_{j} \frac{\lambda}{\lambda_{j}+\lambda}\langle g,v_{j}\rangle v_{j}\,.$ (11) $$

It is known that for classification tasks (what language models do per token), the Gauss-Newton Hessian is equivalent to a form of variance of the generated gradients $\nabla\ell(\hat{y}|x)$, where $\hat{y}\sim p(y|x)$, which is referred to as Fisher Information Matrix (FIM). Generally speaking this is different from the empirical FIM $\mathbf{E}_{z_{train}\sim\mathcal{D}_{train}}\nabla\ell(z_{train})\nabla\ell(z
_{train})^{\top}$, however, for low noise distributions the two might be used interchangeably [Martens, ]. Such interpretation can help us to speculate, that the directions $v_{j}$ corresponding to higher eigenvalues $\lambda_{j}$ are more likely to observe in the training gradients, in the sense that $\mathbf{E}\langle g,v_{j}\rangle^{2}$ is higher. We also notice that $\frac{\lambda}{\lambda+\lambda_{j}}\approx 0$ for the top eigenvalues $\lambda_{j}$ that are much larger than the damping parameter $\lambda$. On the contrary, the lower eigenvalues receive higher weight, that is $\frac{\lambda}{\lambda+\lambda_{j}}\rightarrow 1$ as $\lambda_{j}\rightarrow 0$. In this sense, the iHVP works contrary to the traditional Principal Component Analysis, where the idea is to project the vector onto the top eigenvectors of the covariance. Instead, applying inverse Hessian removes the top directions corresponding to $\lambda_{j}\gg\lambda$ and retains the directions corresponding to $\lambda_{j}\ll\lambda$.

**Table 2: Examples of original and paraphrased sentences.**
| original | “The Mona Lisa, painted by Leonardo Da Vinci in the early 16th century, is one of the most famous paintings in the world.” |
| --- | --- |
| paraphrased | “Leonardo Da Vinci’s early 16th-century painting, the Mona Lisa, is widely regarded as one of the most renowned artworks globally.” |

For example, a plausible interpretation of the directions $v_{j}$ would be that the top directions correspond to general language coherence, sentence structure, and keywords, while the directions with small eigenvalues could correspond to more specific, informative content.
We propose the following experiment to encourage such point of view. We consider ten pairs of sentences, one related to some historical or scientific fact, we refer to as *original*, the other is a paraphrased version of the same fact,
referred to as *paraphrased*.
We show one such pair in Table 2 and in Appendix A we give all 10 pairs(^1^11To avoid cherry-picking, all 20 sentences were generated with Claude 3 Opus with a few prompts.).
For each of the sentences, we calculate the gradient of the next word prediction loss $\nabla\ell(z)$ and calculate pairwise dot-influences $\nabla\ell(z)^{\top}\nabla\ell(z^{\prime})$ and Hessian-based influences $\nabla\ell(z)^{\top}(H+\lambda)^{-1}\nabla\ell(z^{\prime})$.
Our goal is to measure the similarities between original sentences, their rewritings, and their made-up derivatives. For this, we propose to measure the similarity by correspondingly normalizing with norms of gradients and self-influence:

$$ $\displaystyle\text{\rm Gradient-similarity}(z,z^{\prime})$ $\displaystyle=\frac{\nabla\ell(z)^{\top}\nabla\ell(z^{\prime})}{\|\nabla\ell(z )\|\|\nabla\ell(z^{\prime})\|},$ $\displaystyle\text{\rm Influence-similarity}(z,z^{\prime})$ $\displaystyle=\frac{\nabla\ell(z)^{\top}(H+\lambda)^{-1}\nabla\ell(z^{\prime}) }{\sqrt{\nabla\ell(z)^{\top}(H+\lambda)^{-1}\nabla\ell(z)}\sqrt{\nabla\ell(z^{ \prime})^{\top}(H+\lambda)^{-1}\nabla\ell(z^{\prime})}}.$ $$

Figure: Figure 3: Similarity between 20 sentences, see complete list in in Appendix, Section A. Left figure shows influence similarity calculated with LiSSA, middle — gradient similarity, right — the difference between the former and the latter. In the rightmost figure the numbers show the mean over each 10x10 square, with standard error in the brackets. We use OPT 1.3B model, with $\lambda=5.0$, $T=1000$ and $\eta=0.003$. We also use batch size of $4$ sequences, each consisting of $512$ tokens.
Refer to caption: extracted/5879605/figures/simfig_2.png

We show the pairwise similarities between all 20 sentences in Figure 3. In the rightmost graph we also show the difference between gradient similarity and influence similarity. We observe that unrelated sentences generally have higher gradient similarity than influence similarity since the values in the rightmost graph are mostly positive.
As a result, the influence similarity between an original sentence and a rewritten one appears to be consistently higher than between unrelated sentences.

Downweighting directions that are more likely to observe in (11) can also be compared to the idea of the TF-IDF index, where the terms are reweighted according to their inverse frequency [Salton and McGill, ]. Incidentally, we show that for a bag-of-words model (which although trivial, is also a language model), the influence functions correspond to a particular form of the TF-IDF index, see Section D in the Appendix.

## 6 Conclusion

We have shown that the hyperparameters of the classical LiSSA algorithm can be chosen based on two spectral statistics of the Gauss-Newton Hessian: its trace and largest eigenvalue. In particular, we show the batch size used for sampling Hessian-vector products per update has to be sufficiently large. Otherwise, LiSSA might not converge which we demonstrate empirically and theoretically. This particular aspect of hyperparameter choice for LiSSA algorithm has not been previously addressed in the literature. Furthermore, we empirically demonstrate that applying to large models in its original form can still be feasible if we choose a sufficiently large damping parameter.

We do not necessarily advocate for using LiSSA as the algorithm of choice, rather it can be used as a baseline for validating other more lightweight algorithms, in which case it is important to make sure that the hyperparameters of LiSSA are chosen correctly. We hope that our result and the implementation can further facilitate research in influence functions, as well as in other topics where inverse Hessian-vector products naturally appear [Guo et al., , Schulman et al., , Martens, , Kirkpatrick et al., ].

## Appendix A List of prompts for the experiment in Section 5

- 1.
a
original
“The Great Wall of China is the longest wall in the world, stretching over 21,000 kilometers.”
rewrite
“Spanning over 21,000 kilometers, the Great Wall of China holds the record for being the longest wall worldwide.”
- 2.
a
original
“In 1969, Neil Armstrong became the first human to set foot on the Moon during the Apollo 11 mission.”
rewrite
“During the Apollo 11 mission in 1969, Neil Armstrong made history by becoming the first person to walk on the lunar surface.”
- 3.
a
original
“The theory of evolution by natural selection was first proposed by Charles Darwin in his book "On the Origin of Species" in 1859.”
rewrite
“Charles Darwin introduced the concept of evolution by natural selection in his 1859 publication titled "On the Origin of Species".”
- 4.
a
original
“The United Nations was founded in 1945 after World War II to maintain international peace and security.”
rewrite
“Following the conclusion of World War II, the United Nations was established in 1945 to foster international peace and security.”
- 5.
a
original
“The Eiffel Tower, constructed in 1889 for the World’s Fair, is one of the most iconic landmarks in Paris, France.”
rewrite
“One of the most recognizable structures in Paris, France, the Eiffel Tower was built in 1889 for the World’s Fair.”
- 6.
a
original
“The French Revolution, which began in 1789, marked the end of the monarchy and the establishment of a republic in France.”
rewrite
“The monarchy in France was abolished, and a republic was established as a result of the French Revolution, which commenced in 1789.”
- 7.
a
original
“The human brain contains approximately 86 billion neurons, making it the most complex organ in the human body.”
rewrite
“The most intricate organ in the human body, the brain, is composed of roughly 86 billion neurons.”
- 8.
a
original
“In the 2020 United States presidential election, Joe Biden defeated incumbent Donald Trump to become the 46th president.”
rewrite
“Joe Biden secured victory over the sitting president, Donald Trump, in the 2020 United States presidential election, becoming the 46th president.”
- 9.
a
original
“The Mona Lisa, painted by Leonardo da Vinci in the early 16th century, is one of the most famous paintings in the world.”
rewrite
“Leonardo da Vinci’s early 16th-century painting, the Mona Lisa, is widely regarded as one of the most renowned artworks globally.”
- 10.
a
original
“Climate change is a global issue caused by the increase of greenhouse gases in the atmosphere, primarily due to human activities.”
rewrite
“The primary cause of climate change, a worldwide problem, is the accumulation of greenhouse gases in the atmosphere, largely attributed to human activities.”

## Appendix B Postponed proofs from Section 3

### B.1 Proof of Theorem 1

We write,

$$ $\displaystyle\mathbf{u}^{t}-\mathbf{u}^{\star}$ $\displaystyle=\left((1-\lambda\eta)I+\eta\tilde{H}_{t}\right)\mathbf{u}^{t-1}+ \eta\mathbf{g}$ $\displaystyle=\left(I-\eta(\tilde{H}_{t}+\lambda)\right)(\mathbf{u}^{t}-(H+ \lambda)^{-1}\mathbf{g})+\eta(H-\tilde{H}_{t})(H+\lambda)^{-1}\mathbf{g}$ $\displaystyle=\left(I-\eta(\tilde{H}_{t}+\lambda)\right)(\mathbf{u}^{t-1}- \mathbf{u}^{\star})+\eta(H-\tilde{H}_{t})\mathbf{u}^{\star}\,.$ $$

First, taking the expectation and using the fact that $\tilde{H}_{t}$ and $\mathbf{u}^{t}$ are independent,

$$ $\displaystyle\mathbf{E}\mathbf{u}^{t}-\mathbf{u}^{\star}$ $\displaystyle=(I-\eta(H+\lambda))(\mathbf{E}\mathbf{u}^{t-1}-\mathbf{u}^{\star})$ $\displaystyle=(I-\eta(H+\lambda))^{t}(\mathbf{u}^{0}-\mathbf{u}^{\star}),$ $$

and under $\eta<1/\lambda_{\max}(H)$ the matrix $I-\eta(H+\lambda)\preceq(1-\lambda\eta)I$ is a contraction, thus the first bound follows.

For the second part, denote $R=\mathbf{E}(I-\eta(\tilde{H}_{t}+\lambda))^{2}$. Let us take the conditional expectation of the square norm, conditional on all the sampling before step $t$. Setting $\mathcal{F}_{t-1}=\sigma(\tilde{H}_{1},\dots,\tilde{H}_{t-1})$, we have that

$$ $\displaystyle\mathbf{E}\left[\|\mathbf{u}^{t}-\mathbf{u}^{\star}\|^{2}| \mathcal{F}_{t-1}\right]$ $\displaystyle=(\mathbf{u}^{t-1}-\mathbf{u}^{\star})^{\top}\left[\mathbf{E}(I- \eta(\tilde{H}_{t}+\lambda))^{2}\right](\mathbf{u}^{t-1}-\mathbf{u}^{\star})$ $\displaystyle\phantom{=|}+2\eta(\mathbf{u}^{t-1}-\mathbf{u}^{\star})^{\top} \mathbf{E}(I-\eta(\tilde{H}_{t}+\lambda))(H-\tilde{H}^{t})\mathbf{u}^{\star}$ $\displaystyle\phantom{=|}+\eta^{2}\mathbf{E}\|(H-\tilde{H}_{t})\mathbf{u}^{ \star}\|^{2}$ $\displaystyle=\|R^{1/2}(\mathbf{u}^{t-1}-\mathbf{u}^{\star})\|^{2}+2\eta^{2}( \mathbf{u}^{t-1}-\mathbf{u}^{*})^{\top}\{\mathbf{E}\tilde{H}_{t}^{2}-H^{2}\} \mathbf{u}^{\star}$ $\displaystyle\phantom{=|}+\eta^{2}[\mathbf{u}^{\star}]^{\top}\{\mathbf{E} \tilde{H}_{t}^{2}-H^{2}\}\mathbf{u}^{\star}$ $\displaystyle\leq(1-\delta)\|\mathbf{u}^{t-1}-\mathbf{u}^{\star}\|^{2}+2\eta( \mathbf{u}^{t-1}-\mathbf{u}^{\star})^{\top}\{\mathbf{E}\tilde{H}_{t}^{2}-H^{2} \}\mathbf{u}^{\star}+\eta^{2}\tilde{\Delta}\,.$ $$

Here, we have used the fact that by our assumption $\mathbf{E}(I-\eta(\tilde{H_{t}}+\lambda))^{2}\preceq(1-\delta)I$ is a contraction, since we have that

$$ $\displaystyle\mathbf{E}(I-\eta(\tilde{H_{t}}+\lambda))^{2}$ $\displaystyle=(1-\lambda\eta)^{2}I-2\eta(1-\lambda\eta)\mathbf{E}\tilde{H}_{t} +\eta^{2}\mathbf{E}\tilde{H}_{t}^{2}$ $\displaystyle=(1-\lambda\eta)^{2}I-2\eta(1-\lambda\eta){H}+\eta^{2}{H}^{2}$ $\displaystyle\phantom{=|}+\eta^{2}(\mathbf{E}\tilde{H}_{t}^{2}-H^{2})$ $\displaystyle=(1-\eta(H+\lambda))^{2}+\eta^{2}(\mathbf{E}\tilde{H}_{t}^{2}-H^{ 2})\,.$ $$

Taking the unconditional expectation, we obtain that

$$ $\displaystyle\mathbf{E}\|\mathbf{u}^{t}-\mathbf{u}^{\star}\|^{2}$ $\displaystyle\leq(1-\delta)\mathbf{E}\|\mathbf{u}^{t-1}-\mathbf{u}^{\star}\|^{ 2}+2\eta^{2}(\mathbf{u}^{0}-\mathbf{u}^{*})^{\top}(I-\eta(H+\lambda))^{t-1}\{ \mathbf{E}\tilde{H}_{t}^{2}-H\}\mathbf{u}^{\star}+\eta^{2}\tilde{\Delta}$ $\displaystyle\leq\dots$ $\displaystyle\leq(1-\delta)^{t}\|\mathbf{u}^{0}-\mathbf{u}^{\star}\|^{2}+\eta^ {2}\tilde{\Delta}\left(1+(1-\delta)+\dots+(1-\delta)^{t-1}\right)$ $\displaystyle\phantom{=}+2\eta^{2}(\mathbf{u}^{0}-\mathbf{u}^{\star})^{\top} \left\{\sum_{k=0}^{t-1}(1-\delta)^{k}(I-\eta(H+\lambda))^{t-k}\right\}\{ \mathbf{E}\tilde{H}_{t}^{2}-H^{2}\}\mathbf{u}^{\star}$ $\displaystyle\leq(1-\delta)^{t}\|\mathbf{u}^{0}-\mathbf{u}^{\star}\|^{2}+\frac {\eta^{2}\tilde{\Delta}}{1-(1-\delta)}$ $\displaystyle\phantom{=}+2\eta^{2}(\mathbf{u}^{0}-\mathbf{u}^{\star})^{\top} \left\{(1-\delta)^{t}I-(I-\eta(H+\lambda))^{t}\right\}(\eta(H+\lambda)-\delta I )^{-1}\{\mathbf{E}\tilde{H}_{t}^{2}-H^{2}\}\mathbf{u}^{\star}\,.$ $$

We apply the Cauchy-Schwartz inequality to the last term. By (9), we have that $\eta^{2}\{\mathbf{E}\tilde{H}_{t}^{2}-H^{2}\}\preceq\eta(H+\lambda)-\delta I$.
Thus,

$$ $\left\|\left\{(1-\delta)^{t}I-(I-\eta(H+\lambda))^{t}\right\}^{1/2}(\eta(H+ \lambda)-\delta I)^{-1}\eta^{2}\{\mathbf{E}\tilde{H}_{t}^{2}-H^{2}\}\mathbf{u} ^{\star}\right\|^{2}\\ \leq(1-\delta)^{t}\left\|(\eta(H+\lambda)-\delta I)^{-1}\eta^{2}\{\mathbf{E} \tilde{H}_{t}^{2}-H^{2}\}\mathbf{u}^{\star}\right\|^{2}\leq(1-\delta)^{t}\| \mathbf{u}^{\star}\|^{2}.$ $$

We also have,

$$ $\left\|\left\{(1-\delta)^{t}I-(I-\eta(H+\lambda))^{t}\right\}^{1/2}(\mathbf{u} ^{0}-\mathbf{u}^{\star})\right\|^{2}\leq(1-\delta)^{t}\|\mathbf{u}^{0}-\mathbf {u}^{\star}\|^{2}.$ $$

Collecting everything together,

$$ $\mathbf{E}\|\mathbf{u}^{t}-\mathbf{u}^{\star}\|^{2}\leq(1-\delta)^{2}(\| \mathbf{u}-\mathbf{u}^{\star}\|^{2}+\|\mathbf{u}-\mathbf{u}^{\star}\|\|\mathbf {u}^{\star}\|)+\delta^{-1}\eta^{2}\tilde{\Delta}\,.$ $$

### B.2 Proof of Corollary 1

We want to show that with $\frac{C\mathrm{Tr}(H)}{|B|}\geq\eta^{-1}-\lambda$, equation (9) takes place with $\delta=2\eta\lambda-(\eta\lambda)^{2}$. Denote $K=\frac{C\mathrm{Tr}(H)}{|B|}$. Since Condition C.1 holds, we need to show

$$ $(I-\lambda(\eta+H))^{2}+K\eta^{2}H\preceq(1-\delta)I$ $$

The LHS of the above display has eigenvalues $(1-\lambda(\eta+\lambda_{j}))^{2}+K\eta^{2}\lambda_{j}$, where $\lambda_{j}$ are the eigenvalues of $H$. It is therefore sufficient to show that

$$ $\max_{a\in[0,\lambda_{\max}]}(1-\eta(\lambda+a))^{2}+K\eta^{2}a\leq 1-\delta$ $$

We rewrite this condition as

$$ $\max_{a\in[\lambda,\lambda+\lambda_{\max}]}-2\eta a+\eta^{2}a^{2}+K\eta^{2}(a- \lambda)\leq-\delta$ $$

The minimum of the quadratic function is attained at $\bar{a}=\eta^{-1}-K/(2\eta)$. In the case where $\bar{a}$ is in the right half of the interval $[\lambda,\lambda+\lambda_{\max}]$, the maximum of the quadratic function is attained at the point $a=\lambda$. This condition rewrites as $\eta^{-1}-K/(2\eta)\geq\lambda+\lambda_{\max}/2$, which using $\eta(\lambda+\lambda_{\max})=1$ translates into $K\geq\eta^{-1}-\lambda$. Thus, under the assumption that the batch size is at least $|B|\geq\frac{C\mathrm{Tr}(H)}{\eta^{-1}-\lambda}$, we have that

$$ $\max_{a\in[0,\lambda_{\max}]}(1-\eta(\lambda+a))^{2}+K\eta^{2}a=(1-\eta\lambda )^{2}\leq 1-\delta,\qquad\delta=2\eta\lambda-(\eta\lambda)^{2}.$ $$

With such $\lambda$ it holds $(1-\delta)^{t}=(1-\eta\lambda)^{2t}$. It is left to notice that

$$ $\displaystyle\eta^{2}\|\mathbf{E}(\tilde{H}_{t}-H)\mathbf{u}^{\star}\|^{2}$ $\displaystyle=\eta^{2}[\mathbf{u}^{\star}]^{\top}\{\mathbf{E}\tilde{H}_{t}^{2} -H^{2}\}\mathbf{u}^{\star}$ $\displaystyle\leq\frac{(1+c)\mathrm{Tr}(H)}{|B|}\mathbf{g}^{\top}(H+\lambda)^{ -1}H(H+\lambda)^{-1}\mathbf{g}$ $\displaystyle\leq\frac{(1+c)\mathrm{Tr}(H)}{|B|}\mathbf{g}^{\top}(H+\lambda)^{ -1}\mathbf{g}\,.$ $$

### B.3 Motivation for Condition C.1

#### Case of classification.

We first consider the case where the observations in each batch are drawn independently, which could either be the case of classification problem (e.g. on images) or the case where we sample the tokens independently from the set of all training tokens (of course, this would be inefficient and we only consider this case for sake of demonstration).

Consider the gradients

$$ $\mathbf{g}=\nabla\ell(\hat{y}|\mathbf{x}),\qquad\mathbf{x}\sim\mathcal{D}_{tr} ,\;\;\hat{y}|\mathbf{x}\sim p(\mathbf{x}).$ $$

Here, $\mathbf{x}\sim\mathcal{D}_{tr}$ either means that $\mathbf{x}$ is drawn uniformly from a finite training dataset, or from some population distribution. The expression $\hat{y}|\mathbf{x}\sim p(\mathbf{x})$ means that once we have drawn $\mathbf{x}$, we then generate a random label $\hat{y}$ according to the model’s output distribution $p(\mathbf{x})$. For language models, $\mathbf{x}=(s_{1},\dots,s_{t-1})$ is the context and $\hat{y}$ predicts $s_{t}$. For simplicity, we assume that $\mathbf{g}$ are spherically distributed, meaning that there are independent $r$, $\mathbf{s}$ such that $\mathbf{E}r^{2}=1$ and $\mathbf{s}$ is distributed uniformly over a unit sphere in $\mathbb{R}^{N}$, and $\mathbf{g}\overset{d}{=}\sqrt{N}rH^{1/2}\mathbf{s}$. This helps to account for heavy tails of the gradients that are often observed in practice, which would not be captured, for example, by sub-Gaussian distributions.

The in-batch Gauss-Newton Hessian then reads as

$$ $\tilde{H}_{t}=\frac{1}{|B|}\sum_{\mathbf{x}\in B}\mathbf{E}_{\hat{y}\sim p( \mathbf{x})}\nabla\ell(\hat{y}|\mathbf{x})\nabla\ell(\hat{y}|\mathbf{x})^{\top}$ $$

Thanks to the fact that the elements in $B$ are i.i.d. we have that

$$ $\displaystyle\mathbf{E}\tilde{H}_{t}^{2}=\,$ $\displaystyle\frac{1}{|B|}\mathbf{E}\sum_{\mathbf{x}\neq\mathbf{x}^{\prime}} \mathbf{E}_{\hat{y}\sim p(\mathbf{x})}\nabla\ell(\hat{y}|\mathbf{x})\nabla\ell (\hat{y}|\mathbf{x})^{\top}\mathbf{E}_{\hat{y}\sim p(\mathbf{x}^{\prime})} \nabla\ell(\hat{y}|\mathbf{x}^{\prime})\nabla\ell(\hat{y}|\mathbf{x}^{\prime}) ^{\top}$ $\displaystyle+\frac{1}{|B|^{2}}\sum_{\mathbf{x}\in B}\left[\mathbf{E}_{\hat{y} \sim p(\mathbf{x})}\nabla\ell(\hat{y}|\mathbf{x})\nabla\ell(\hat{y}|\mathbf{x} )^{\top}\right]^{2}$ $\displaystyle=\,$ $\displaystyle(1-1/|B|)H^{2}+\frac{1}{|B|^{2}}\sum_{\mathbf{x}\in B}\left[ \mathbf{E}_{\hat{y}\sim p(\mathbf{x})}\nabla\ell(\hat{y}|\mathbf{x})\nabla\ell (\hat{y}|\mathbf{x})^{\top}\right]^{2}$ $\displaystyle\preceq\,$ $\displaystyle(1-1/|B|)H^{2}+\frac{1}{|B|}\mathbf{E}\|\nabla\ell(\hat{y}| \mathbf{x})\|^{2}\nabla\ell(\hat{y}|\mathbf{x})\nabla\ell(\hat{y}|\mathbf{x})^ {\top}$ $$

Now, writing $\mathbf{g}$ instead of $\nabla\ell(\hat{y}|\mathbf{x})$ and using the spherical property $\mathbf{g}=\sqrt{N}rH^{1/2}\mathbf{s}$, we have that

$$ $v^{\top}[\mathbf{E}\|\mathbf{g}\|^{2}\mathbf{g}\mathbf{g}^{\top}]v=N^{2} \mathbf{E}r^{4}\mathbf{E}\mathbf{s}^{\top}H\mathbf{s}(\mathbf{s}^{\top}H^{1/2} vv^{\top}H^{1/2}\mathbf{s})$ $$

For $\mathbf{s}$ uniformly distributed over a unit sphere, we can use the identity $\mathbf{E}(\mathbf{s}^{\top}A\mathbf{s})(\mathbf{s}^{\top}B\mathbf{s})=\frac{1
}{N(N+2)}(2\mathrm{Tr}(AB)+\mathrm{Tr}(A)\mathrm{Tr}(B))$ [Wiens, ]. Then, we have that

$$ $v^{\top}[\mathbf{E}\|\mathbf{g}\|^{2}\mathbf{g}\mathbf{g}^{\top}]v=(\mathbf{E} r^{4})\frac{N}{N+2}(2v^{\top}H^{2}v+\mathrm{Tr}(H)v^{\top}Hv),$ $$

therefore,

$$ $\mathbf{E}\tilde{H}_{t}^{2}-H^{2}\preceq\frac{1}{|B|}\left\{(\mathbf{E}r^{4}) \mathrm{Tr}(H)H+(\mathbf{E}r^{4}-1)H^{2}\right\}.$ $$

#### Case of incoherent tokens.

Let us now consider the case where we sample the tokens sequence-wise. For simplicity, we assume that each sequence has the same size, so that each token in the dataset has equal probability to be drawn (when we say token to be drawn, we mean that we consider the prediction of token $s_{t}$ with context $(s_{1},\dots,s_{t-1})$. This means that when $s\in B$, and $|s|=L$, our batch contains all contexts $(s_{1},\dots,s_{t-1})$ for $t=1,\dots,L$. Here, we assume that gradients corresponding to tokens in the same sequence have little correlations, a condition that following Tang et al. [] we call *incoherence*. The idea is that for high-dimensional problems it is natural to expect that gradients sampled in the same batch are likely to have low cosine similarity.

For simplicity, we formulate it pointwise for every sequence $s$ in the population. Assume that there is a bounded $R\geq 1$, such that for all $s\in\mathcal{D}_{tr}$, for all $\mathbf{x}\in s$,

$$ $\displaystyle\sum_{\mathbf{x}^{\prime}\in s}\max_{\hat{y},\hat{y}^{\prime}}| \cos(\nabla\ell(\hat{y},\mathbf{x}),\nabla\ell(\hat{y}^{\prime},\mathbf{x}^{ \prime}))|\leq R\,.$ $$

One way to interpret such condition is as follows: if we look at the cosine similarity as a measure of relevance between the tokens that we predict, we can say that in each sequence there is a limited number of keywords which relate to each other, while the rest of the words are only related locally.

Similar to the independent case, we also assume that for a randomly sampled token, $\nabla\ell(\hat{y}|\mathbf{x})$ has a spherical distribution, i.e. $\nabla\ell(\hat{y}|\mathbf{x})\overset{d}{=}\sqrt{N}rH^{1/2}\mathbf{s}$.

Let $B=s_{1}\cup\dots\cup s_{b}$, where $b$ is the number of sequences in a batch, so that $|B|=bL$, and we assume that each sequence has the same length $L+1$ (we do not predict the first token in a sequence, whose index is 0 0). Let $\tilde{H}(s)$ denotes in-sequence GNH. Then,

$$ $\mathbf{E}\tilde{H}_{t}^{2}=H^{2}+\frac{1}{b}\{\mathbf{E}\tilde{H}(s)^{2}-H^{2 }\},$ (12) $$

where $s$ is a single random sequence. Let us expand,

$$ $\tilde{H}(s)^{2}=\frac{1}{L^{2}}\sum_{\mathbf{x},\mathbf{x}^{\prime}\in s} \tilde{H}(\mathbf{x})\tilde{H}(\mathbf{x}^{\prime})$ $$

We have that for $\mathbf{x},\mathbf{x}^{\prime}\in s$,

$$ $\displaystyle\frac{1}{2}$ $\displaystyle(\tilde{H}(\mathbf{x})\tilde{H}(\mathbf{x}^{\prime})+\tilde{H}( \mathbf{x}^{\prime})\tilde{H}(\mathbf{x}))$ $\displaystyle=\frac{1}{2}\mathbf{E}_{\hat{y}\sim p(\mathbf{x}),\hat{y}^{\prime }\sim p(\mathbf{x}^{\prime})}\langle\nabla\ell(\hat{y}|\mathbf{x}),\nabla\ell( \hat{y}^{\prime}|\mathbf{x}^{\prime})\rangle(\nabla\ell(\hat{y}|\mathbf{x}) \nabla\ell(\hat{y}^{\prime}|\mathbf{x}^{\prime})^{\top}+\nabla\ell(\hat{y}^{ \prime}|\mathbf{x}^{\prime})\nabla\ell(\hat{y}|\mathbf{x})^{\top})$ $\displaystyle\preceq\frac{\delta(\mathbf{x},\mathbf{x}^{\prime})}{2}\left( \mathbf{E}_{\hat{y}\sim p(\mathbf{x})}\|\nabla\ell(\hat{y}|\mathbf{x})\|^{2} \nabla\ell(\hat{y}|\mathbf{x})\nabla\ell(\hat{y}|\mathbf{x})^{\top}+\mathbf{E} _{\hat{y}^{\prime}\sim p(\mathbf{x}^{\prime})}\|\nabla\ell(\hat{y}^{\prime}| \mathbf{x}^{\prime})\|^{2}\nabla\ell(\hat{y}^{\prime}|\mathbf{x}^{\prime}) \nabla\ell(\hat{y}^{\prime}|\mathbf{x}^{\prime})^{\top}\right),$ $$

where we denote for short $\delta(\mathbf{x},\mathbf{x}^{\prime})=\max_{\hat{y},\hat{y}^{\prime}}|\cos(
\nabla\ell(\hat{y},\mathbf{x}),\nabla\ell(\hat{y}^{\prime},\mathbf{x}^{\prime}
))|$, and we also
use the fact that $2\|a\|\|b\|\langle v,a\rangle\langle v,b\rangle\leq\|a\|^{2}\langle v,a\rangle
^{2}+\|b\|^{2}\langle v,b\rangle^{2}$.
Summing up we have that

$$ $\displaystyle\frac{1}{L^{2}}\sum_{\mathbf{x},\mathbf{x}^{\prime}\in s}\tilde{H }(\mathbf{x})\tilde{H}(\mathbf{x}^{\prime})$ $\displaystyle\preceq\frac{1}{L^{2}}\sum_{\mathbf{x}\in s}\left(\sum_{\mathbf{x }^{\prime}\in s}\delta(\mathbf{x},\mathbf{x}^{\prime})\right)\mathbf{E}_{\hat{ y}\sim p(\mathbf{x})}\|\nabla\ell(\hat{y}|\mathbf{x})\|^{2}\nabla\ell(\hat{y}| \mathbf{x})\nabla\ell(\hat{y}|\mathbf{x})^{\top}$ $\displaystyle\preceq\frac{R}{L^{2}}\sum_{\mathbf{x}\in s}\mathbf{E}_{\hat{y} \sim p(\mathbf{x})}\|\nabla\ell(\hat{y}|\mathbf{x})\|^{2}\nabla\ell(\hat{y}| \mathbf{x})\nabla\ell(\hat{y}|\mathbf{x})^{\top}$ $$

From previous section we know that for $\mathbf{g}=\nabla\ell(\hat{y}|\mathbf{x})$, we have $\mathbf{E}\|\mathbf{g}\|^{2}\mathbf{g}\mathbf{g}^{\top}\preceq(\mathbf{E}r^{4}
)(\mathrm{Tr}(H)H+H^{2})$. Taking the expectation over $s\sim\mathcal{D}_{tr}$ and $\mathbf{x}\sim\mathcal{U}(s)$ we obtain,

$$ $\mathbf{E}\tilde{H}(s)^{2}\preceq\frac{R}{L}(\mathbf{E}r^{4})(\mathrm{Tr}(H)H+ H^{2})$ $$

Plugging into (12), we obtain that

$$ $\mathbf{E}\tilde{H}_{t}^{2}-H^{2}\preceq\frac{R}{|B|}(\mathbf{E}r^{4})\mathrm{ Tr}(H)H+\left(\frac{R}{|B|}(\mathbf{E}r^{4})-\frac{1}{b}\right)H^{2}\,.$ $$

#### Empirical sanity check.

Furthermore we propose a simple empirical sanity check by assessing the relationship between the LHS and RHS in condition (C.1). In Figure 4, we compare the traces of two matrices, similarly evaluating the traces by averaging over random quadratic forms $\mathbf{g}^{\top}(\tilde{H}_{t}^{2}-H^{2})\mathbf{g}$. We evaluate this gap by estimating the HVPs $\tilde{H}_{t}\mathbf{g}$, on random Gaussian vectors and taking their norm. Each evaluation is averaged over 1000 realizations of $\mathbf{g}$.

Figure: Figure 4: Compring traces of LHS and RHS of condition C.1 for different batch sizes. We evaluate the traces for ResNet-18, ResNet-50, and OPT-1.3B, 4 batch sizes for each model. For the OPT-1.3B, the batch size is counted in tokens.
Refer to caption: extracted/5879605/figures/condition1.png

### B.4 Counter-example with divergence

###### Lemma 1 .

There exists a binary classification task, where condition C.1 holds with exact equality

$$ $\mathbf{E}\tilde{H}_{t}^{2}-H^{2}=\frac{1}{|B|}\mathrm{Tr}(H)H,$ $$

and for some inputs $\mathbf{u}^{0},\mathbf{g}$ the following claim holds.

Let us choose arbitrary step size satisfying $0<\eta\leq 1/\lambda_{\max}(H)$, such the matrix

$$ $(1-\eta(H+\lambda))^{2}+\eta^{2}\{\mathbf{E}\tilde{H}_{t}^{2}-H^{2}\}$ $$

has eigenvalue $\lambda\geq 1$. Then, the LiSSA algorithm converges on average, but not samplewise:

$$ $\|\mathbf{E}\mathbf{u}^{t}-\mathbf{u}^{\star}\|\rightarrow 0,\qquad\mathbf{E} \|\mathbf{u}^{t}-\mathbf{u}^{\star}\|^{2}\rightarrow\infty\,.$ $$

In particular, such $\eta$ exists in the case $|B|\leq\left(\frac{\lambda_{\max}}{\lambda_{\max}+\lambda}\right)^{2}\frac{
\mathrm{Tr}(H)}{\lambda_{\max}}$.

###### Proof.

We consider a binary regression with $y\in\{0,1\}$ and inputs $\mathbf{x}\in\mathbb{R}^{N}$. For simplicity, we assume that $\mathbf{x},y$ are independent, and $y$ takes values $\{0,1\}$ with equal probabilities. We consider the distribution in $\mathbf{x}$ that satisfies $\|\mathbf{x}\|^{2}=\mathrm{Tr}(H)$ pointwise and $\mathbf{E}\mathbf{x}\mathbf{x}^{\top}=H$.
For this, let $H=V\Lambda V^{\top}$ where $\Lambda=\mathrm{diag}\{\lambda_{1},\dots,\lambda_{N}\}$ and take $\mathbf{x}=V\mathbf{s}$, $\mathbf{s}=(\sqrt{\lambda_{1}}\epsilon_{1},\dots,\sqrt{\lambda_{N}}\epsilon_{N
})^{\top}$ where $\epsilon_{i}=\pm 1$ with equal probabilities. Consider the binary logistic model

$$ $\displaystyle\log p(y=1|\mathbf{x})$ $\displaystyle=\theta^{\top}\mathbf{x}-\log(\exp(\theta^{\top}\mathbf{x})+\exp( -\theta^{\top}\mathbf{x})),$ $\displaystyle\log p(y=0|\mathbf{x})$ $\displaystyle=-\theta^{\top}\mathbf{x}-\log(\exp(\theta^{\top}\mathbf{x})+\exp (-\theta^{\top}\mathbf{x})),$ $$

and assume that during training the model converged to the optimal parameters $\theta=0$, since $\mathbf{x},y$ are independent. Then it is straightforward to calculate that

$$ $\displaystyle\tilde{H}_{t}=\frac{1}{|B|}\sum_{(\mathbf{x},y)\in B}\mathbf{x} \mathbf{x}^{\top},\qquad\mathbf{E}\tilde{H}_{t}^{2}=\left(1-\frac{1}{|B|} \right)H^{2}+\frac{1}{|B|}\mathrm{Tr}(H)H.$ $$

Since all matrices are rotated by $V$ from left and right, we assume w.l.o.g. that $V=I$.

Now assume that $\mathbf{g}=0$ we have $\mathbf{u}^{\star}=0$, and let $\mathbf{u}^{0}\neq 0$. We have

$$ $\mathbf{u}^{t}=(1-\eta(\tilde{H}_{t}+\lambda))\mathbf{u}^{t-1}=\prod_{j=1}^{t} (1-\eta(\tilde{H}_{j}+\lambda))\mathbf{u}^{0}\,.$ $$

Set $Q=I-\eta(\hat{H}_{t}+\lambda)$ and $R=\mathbf{E}(1-\eta(\hat{H}_{t}+\lambda))^{2}=(1-\eta(H+\lambda))^{2}+\eta^{2}
(\mathbf{E}\tilde{H}_{t}^{2}-H^{2})$. Notice that both matrices are diagonal. Consider the sequence of matrices $R_{0}=I$, $R_{1}=R$, $R_{k}=\mathbf{E}(1-\eta(\tilde{H}_{t}+\lambda))R_{k-1}(1-\eta(\tilde{H}_{t}+
\lambda))$. Also denote partial product $\tilde{T}_{j}=\prod_{k=1}^{j}(1-\eta(\tilde{H}_{k}+\lambda))$. Then we have

$$ $\displaystyle\mathbf{E}\|\mathbf{u}^{t}\|^{2}$ $\displaystyle=[\mathbf{u}^{0}]^{\top}\mathbf{E}C_{t-1}(1-\eta(\tilde{H}_{t}+ \lambda))^{2}C_{t-1}\mathbf{u}^{0}$ $\displaystyle=[\mathbf{u}^{0}]^{\top}\mathbf{E}T_{t-1}RC_{t-1}\mathbf{u}^{0}$ $\displaystyle=[\mathbf{u}^{0}]^{\top}\mathbf{E}T_{t-2}(1-\eta(\tilde{H}_{t-1}+ \lambda))R_{1}(1-\eta(\tilde{H}_{t-2}+\lambda))C_{t-2}\mathbf{u}^{0}$ $\displaystyle=[\mathbf{u}^{0}]^{\top}\mathbf{E}T_{t-2}R_{2}T_{t-2}\mathbf{u}^{0}$ $\displaystyle=[\mathbf{u}^{0}]^{\top}\mathbf{E}T_{t-3}R_{3}T_{t-3}\mathbf{u}^{0}$ $\displaystyle=\dots$ $\displaystyle=[\mathbf{u}^{0}]^{\top}R_{t}\mathbf{u}^{0}$ $$

Let us show that the matrix $B_{t}$ is diagonal. Indeed, assuming $B_{k-1}$ is diagonal, we have that

$$ $\displaystyle R_{k}$ $\displaystyle=QR_{k-1}Q+\eta^{2}\mathbf{E}(\tilde{H}_{t}-H)R_{k-1}(\tilde{H}_{ t}-H)$ $\displaystyle=QB_{k-1}Q-\eta^{2}HR_{k-1}H+\eta^{2}\mathbf{E}\tilde{H}_{t}R_{k- 1}\tilde{H}_{t}$ $$

We have that for diagonal $R_{k-1}$, $\mathbf{x}^{\top}R_{k-1}\mathbf{x}=\sum_{j}\lambda_{j}R_{k-1}[j,j]=Tr(HR_{k-1})$ is deterministic. Therefore,

$$ $\displaystyle\mathbf{E}\left(\frac{1}{|B|}\sum_{x\in B}\mathbf{x}\mathbf{x}^{ \top}\right)R_{k-1}\left(\frac{1}{|B|}\sum_{x\in B}\mathbf{x}\mathbf{x}^{\top}\right)$ $\displaystyle=\left(1-\frac{1}{|B|}\right)H^{2}R_{k-1}+\frac{1}{|B|}\mathbf{E} (\mathbf{x}^{\top}R_{k-1}\mathbf{x})\mathbf{x}\mathbf{x}^{\top}$ $\displaystyle=\left(1-\frac{1}{|B|}\right)H^{2}R_{k-1}+\mathrm{Tr}(HR_{k-1})H$ $\displaystyle=[\mathbf{E}\tilde{H}_{t-1}^{2}]R_{k-1}$ $$

Thus, we conclude

$$ $R_{k}=Q^{2}R_{k-1}-\eta^{2}H^{2}R_{k-1}+\eta^{2}[\mathbf{E}\tilde{H}_{t}^{2}]R _{k-1}=RR_{k-1}=R^{k}R_{0}=R^{k}.$ $$

Now, we have that

$$ $\mathbf{E}\|\mathbf{u}^{t}\|^{2}=[\mathbf{u}^{0}]^{\top}R^{t}\mathbf{u}^{0} \geq\lambda^{t}\langle\mathbf{u}^{0},v\rangle^{2}\not\to 0,$ $$

whenever $\langle\mathbf{u}^{0},v\rangle\neq 0$.
∎

###### Lemma 1 .

###### Proof.

## Appendix C Details of Hessian statistics calculation

For trace approximation, we use the fact that for a Gaussian vector $\mathbf{g}\in\mathcal{N}(0,\tfrac{1}{N}I)$, generated independently from $\tilde{H}$, $\mathbf{E}\mathbf{g}^{\top}\tilde{H}\mathbf{g}=\frac{1}{N}\mathrm{Tr}(H)\,$.
Observe that from (7),

$$ $\mathbf{g}^{\top}\tilde{H}\mathbf{g}=\sum_{j=1}^{N}\left\{J_{\theta}h(\mathbf{ x};\theta)^{\top}\mathbf{g}\right\}^{\top}\{\mathrm{Diag}(\mathrm{sf}(h))- \mathrm{sf}(h)\mathrm{sf}(h)^{\top}\}[J_{\theta}h(\mathbf{x};\theta)^{\top} \mathbf{g}],$ $$

thus evaluating such quadratic form can be done by only evaluating the middle Hessian $\mathrm{Diag}(\mathrm{sf}(h))-\mathrm{sf}(h)\mathrm{sf}(h)^{\top}$ and finite differences $J_{\theta}h(\mathbf{x};\theta)^{\top}\mathbf{g}\approx 50(h(\mathbf{x};\theta+
0.01\mathbf{g})-h(\mathbf{x};\theta-0.01\mathbf{g}))$. The results reported in Table 1 are based on 1600 evaluations of $\mathbf{g}_{i}^{\top}\tilde{H}_{i}\mathbf{g}_{i}$, including the standard error.

At an increased price, we can also evaluate the Frobenius norm of the Hessian. For that, for a given vector $\mathbf{g}$, sampled from Gaussian distribution, we evaluate two independent Hessian-vector products $\tilde{H}\mathbf{g}$, $\hat{H}\mathbf{g}$, so that $\frac{1}{N}\mathrm{Tr}(H^{2})=\mathbf{E}(\tilde{H}\mathbf{g})^{\top}\hat{H}
\mathbf{g}$. Then, by sampling independently a series of independent realizations, we can evaluate the mean $\frac{1}{N}\mathrm{Tr}(H^{2})$ and the standard error of our estimation. We report these evaluations in Table 3.

For evaluating the top eigenvalues, one can employ the sketching technique. For example, Swartworth and Woodruff [] show that for $\Phi\in\mathbb{R}^{d\times N}$ generated in a way such that $\Phi_{ij}\sim\mathcal{N}(0,\tfrac{1}{d})$ for $d$ large enough, we have that for finite amount of top eigenvalues,

$$ $\lambda_{l}(H)\approx\lambda_{l}(\Phi H\Phi^{\top}).$ (13) $$

In order to evaluate the matrix $\Phi H\Phi^{\top}$, we iterate over each of $d$ columns ${\boldsymbol{\phi}}_{j}$ of $\Phi^{\top}$ and evaluate the HVP $H{\boldsymbol{\phi}}$ by sampling empirical $\tilde{H}_{t}{\boldsymbol{\phi}}$. For language models, to speed up the caluculation we truncate the context length to $256$ and average over batch of size $50$. For image classification, we use a batch size of $5000$ to evaluate each $H{\boldsymbol{\phi}}$. We then project each of these columns back with the embedding $\Phi$, so the result is a $d\times d$ matrix, whose top eigenvalue can be calculated with standard linear algebra packages. According to Swartworth and Woodruff [],
the error term in (13) is bounded by multiple of $\|H\|_{Fr}/\sqrt{d}$. Although we do not offer a precise control, we suggest that taking $d=5000$ should be sufficient for the models in Table 3. Note that their bound does not account for sampling error and we leave it out of consideration in the scope of this paper, we simply want to produce some adequate bound on the largest eigenvalue.

**Table 3: Second moment statistics for Gauss-Newton-Hessian calculated on ImageNet (IN) on vision and Open-Web-Text-2 (OWT) on language models.**
| Model | Size | Data | $\|H\|_{Fr}$ | $\lambda_{\max}(H)$ | $\left(\frac{\|H\|_{Fr}}{\lambda_{\max}(H)}\right)^{2}$ |
| --- | --- | --- | --- | --- | --- |
| ResNet-18 | 11M | IN | $2.55\times 10^{4}$ | $\approx 270$ | $8.92\times 10^{3}$ |
| ResNet-50 | 25M | IN | $2.67\times 10^{4}$ | $\approx 470$ | $3.22\times 10^{3}$ |
| OPT | 1.3B | OWT | $2.94\times 10^{3}$ | $\approx 780$ | $1.42\times 10^{1}$ |
| Llama-1 | 7B | OWT | $3.73\times 10^{3}$ | $\approx 1600$ | $5.43$ |
| Mistral | 7B | OWT | $2.49\times 10^{4}$ | $\approx 5600$ | $1.98\times 10^{1}$ |

#### Pseudo-random embeddings.

Generating and storing a dense embedding matrix $\Phi\in\mathbb{R}^{d\times N}$ can be prohibitively expensive, since each of the rows $\Phi[i,:]$ is equivalent to one model in memory. We propose instead to use pseudo-random generators, and generate $\Phi$ “on-the-fly” using a single integer number *seed*. The model parameters are usually accessed in a form of lists $\theta^{\top}=(\theta_{1}^{\top},\dots,\theta_{L}^{\top})$ say where $L$ is the number of layers. Correspondingly, the gradients and HVPs are also iterated over a list $\mathbf{g}=(\mathbf{g}_{1},\dots,\mathbf{g}_{L})$. We generate $\Phi$ in the form $\Phi=\begin{pmatrix}\Phi_{1}&\dots&\Phi_{L}\end{pmatrix}$, so that $\Phi\mathbf{g}=\sum_{j=1}^{L}\Phi_{j}\mathbf{g}_{j}$. To calculate $\Phi\mathbf{g}$ we initialize a random generator with a fixed seed, and then generate $\Phi_{j}$ on the fly, so that we never have to store the whole matrix $\Phi$ in the memory.

Nevertheless, the embedding operation itself $\Phi\mathbf{g}$ can be too expensive. Our observation is that starting $d=50$ we do not really benefit from parallel matrix computations on GPU and the price scales linearly, i.e. $d=200$ is $4$ times as expensive as $d=50$, and so on. Furthermore, the price of a single application of embedding for $d=50$ can be as high as the gradient computation itself. In order to reduce the computational price of embedding, we suggest to use the following heuristic.
Instead of summing up per-layer embeddings
$\Phi\mathbf{g}=\Phi_{1}\mathbf{g}_{1}+\dots+\Phi_{L}\mathbf{g}_{L}$, we suggest to concatenate them.
This way, we can increase the dimension by a factor of $L$ with little additional computation. Let us denote the resulting concatenating embedding by ${\boldsymbol{\{\Phi\}}}$, then we can write

$$ $\displaystyle\Phi$ $\displaystyle=\begin{pmatrix}\Phi_{1}&\Phi_{2}&\dots&\Phi_{L}\end{pmatrix}\in \mathbb{R}^{d\times N}$ $\displaystyle{\boldsymbol{\{\Phi\}}}$ $\displaystyle=\begin{pmatrix}\Phi_{1}&0&\dots&\\ 0&\Phi_{1}&\\ &&\dots&\\ 0&0&\dots&\Phi_{L}\end{pmatrix}\in\mathbb{R}^{Ld\times N}$ $$

In other words, we have ${\boldsymbol{\{\Phi\}}}\theta=\mathtt{vstack}([\Phi_{1}\theta_{1},\dots,\Phi_{
L}\theta_{L}])$. That is, each $\Phi_{1}\theta_{1}$ has dimension $d$ and

$$ $\mathbf{E}_{\Phi}[{\boldsymbol{\{\Phi\}}}\theta]^{\top}[{\boldsymbol{\{\Phi\}} }\theta^{\prime}]=\sum_{j}\mathbf{E}[\Phi_{j}\theta_{j}]^{\top}[\Phi_{j}\theta _{j}^{\prime}]=\sum_{j}\theta_{j}^{\top}\theta_{j}^{\prime}=\theta^{\top} \theta^{\prime},$ $$

so that on average it preserves the dot products as well. However, concatenation can dramatically reduce the variance of one dot product. Indeed, we have that

$$ $\displaystyle\mathrm{Var}_{\Phi}([{\boldsymbol{\{\Phi\}}}\theta]^{\top}[{ \boldsymbol{\{\Phi\}}}\theta^{\prime}])$ $\displaystyle=\sum_{j}\mathrm{Var}([\Phi_{j}\theta_{j}]^{\top}[\Phi_{j}\theta_ {j}^{\prime}])\lesssim\frac{1}{d}\sum_{j}\|\theta_{j}\|^{2}\|\theta_{j}^{ \prime}\|^{2},$ $$

and recall the bound from before,

$$ $\mathrm{Var}_{\Phi}([\Phi\theta]^{\top}[\Phi\theta^{\prime}])\lesssim\frac{1}{ d}\|\theta\|^{2}\|\theta^{\prime}\|^{2}=\frac{1}{d}\left(\sum_{j}\|\theta_{j} \|^{2}\right)\left(\sum_{j}\|\theta_{j}^{\prime}\|^{2}\right)$ $$

The former can be much smaller than the latter when square norms of the gradients are spread “evenly” over the layers. That is, assume that $\|\theta_{j}\|^{2}$ is approximately in the same bulk $C^{-1}M\leq\|\theta_{j}\|^{2}\leq CM$, $M=\frac{1}{L}\sum_{j}\|\theta_{j}\|^{2}$. Then $\mathrm{Var}_{\Phi}([{\boldsymbol{\{\Phi\}}}\theta]^{\top}[{\boldsymbol{\{\Phi
\}}}\theta^{\prime}])\lesssim\frac{L}{d}M^{2}$ while $\mathrm{Var}_{\Phi}([\Phi\theta]^{\top}[\Phi\theta^{\prime}])\lesssim\frac{L^{
2}}{d}M^{2}$, so we effectively it is equivalent to increasing dimension $L$ times with original Gaussian features.
Here we ignored the terms $\langle\theta,\theta^{\prime}\rangle^{2}$ but in practice they are significantly smaller than $\|\theta\|^{2}\|\theta^{\prime}\|^{2}$.

## Appendix D Influence functions for bag-of-words model is TF-IDF

Recall that for a set of documents $d\in D$ and terms $t\in T$. Define term frequency (TF) and document frequency (DF) as follows,

$$ $\displaystyle TF(t,d)$ $\displaystyle=\frac{1}{|d|}count(t,d),$ $\displaystyle DF(t)$ $\displaystyle=\frac{\#\left\{d\in D:t\in d\right\}}{|D|};$ $$

let us also consider a variant of inverse document frequency $IDF(t)=\sqrt{1/DF(t)}$. Note that in standard $TF{\cdot}IDF$ definition, the square root is replaced with logarithm, we only propose to consider the square root for sake of comparison to influence functions. Then, the document-term relevance $TF\cdot IDF$ is calculated as the product of $TF$ and $IDF$, and the corresponding similarity between $d_{1}$ and $d_{2}$ is

$$ $\mathrm{sim}(d_{1},d_{2})=\sum_{t}TF\cdot IDF(d_{1},t)TF\cdot IDF(d_{2},t)= \sum_{t}\frac{1}{DF(t)}TF(d_{1},t)TF(d_{2},t)$ $$

The bag-of-words model reads as follows,

$$ $\log p(d)=\sum_{t}\mathtt{count}(t,d)\log p_{t},$ $$

where $p_{t}=\exp(x_{t})/(\exp(x_{1})+\dots+\exp(x_{T}))$, the $x_{t}$ are parameters. For sake of simplicity we assume that each document has the same length $|d|$. We have,

$$ $\displaystyle\nabla_{x}\log p(d)$ $\displaystyle=\sum_{t}{count}(t,d)e_{t}-|d|\nabla_{x}\log\left(\sum_{t}\exp(x_ {t})\right)=\sum_{t}{count}(t,d)e_{t}-|d|\mathrm{sf}(x)$ $\displaystyle=|d|\left\{\sum_{t}TF(t,d)e_{t}-\mathrm{sf}(x)\right\}$ (14) $$

Notice that by definition, $\mathbf{1}^{\top}\nabla_{x}\log p(d)=0$.
The Hessian looks as follows

$$ $\nabla_{x}^{2}\mathbf{E}\log p(d)=-|d|\nabla_{x}^{\top}\mathrm{sf}(x)=-|d| \left(\mathrm{Diag}(\mathrm{sf}(x))-\mathrm{sf}(x)\mathrm{sf}(x)^{\top}\right),$ $$

where we calculate that

$$ $\nabla_{x_{i}}\frac{\exp(x_{j})}{\sum_{k}\exp(x_{k})}=\frac{\delta_{ij}\exp(x_ {i})(\sum_{k}\exp(x_{k}))-\exp(x_{i})\exp(x_{j})}{(\sum_{k}\exp(x_{k}))^{2}}= \delta_{ij}\mathrm{sf}(x)_{i}-\mathrm{sf}(x)_{i}\mathrm{sf}(x)_{j}$ $$

Let us calculate the inverse Hessian. Given a damping parameter $\lambda$, let $\mathbf{q}=(\mathbf{p}+\lambda)^{1/2}$ elementwise. Note that $\|\mathbf{q}\|^{2}=1+\lambda N$. Set also $\mathbf{r}=\mathbf{p}/(\mathbf{p}+\lambda)^{1/2}$ elementwise, and notice that $\|r\|^{2}=\sum_{j}p_{j}^{2}/(p_{j}+\lambda)=1-\lambda\sum_{j}p_{j}/(p_{j}+
\lambda)<1$. Then,

$$ $\mathrm{Diag}(\mathbf{p})+\lambda I-\mathbf{p}\mathbf{p}^{\top}=\mathrm{Diag}( \mathbf{q})(I-\mathbf{r}\mathbf{r}^{\top})\mathrm{Diag}(\mathbf{q}),$ $$

so the inverse equals to

$$ $\displaystyle(H+\lambda)^{-1}$ $\displaystyle=\mathrm{Diag}(\mathbf{q})^{-1}\left(I+\frac{1}{1-\|\mathbf{r}\|^ {2}}\mathbf{r}\mathbf{r}^{\top}\right)\mathrm{Diag}(\mathbf{q})^{-1}$ $\displaystyle=\mathrm{Diag}(\mathbf{p}+\lambda)^{-1}-\left(\lambda\sum_{j}p_{j }/(p_{j}+\lambda)\right)^{-1}\mathbf{d}\mathbf{d}^{\top},$ $$

where we denote $\mathbf{d}=\mathbf{p}/(\mathbf{p}+\lambda)=\mathbf{1}-\lambda/(\mathbf{p}+\lambda)$ elementwise.
Note that the gradients $\mathbf{g}_{d}=\nabla\log p(d)$ are in the subset $\mathbf{g}_{d}^{\top}\mathbf{1}$. Since $\mathbf{d}=\mathbf{1}+O(\lambda)$ we have that

$$ $\mathbf{g}_{d_{1}}^{\top}\left(\lambda\sum_{j}p_{j}/(p_{j}+\lambda)\right)^{-1 }\mathbf{d}\mathbf{d}^{\top}\mathbf{g}_{d_{2}}=O(\lambda),$ $$

therefore, in the limit $\lambda\rightarrow 0$, we have that the influence between documents $d_{1},d_{2}$ reads as

$$ $\mathcal{I}(d_{1},d_{2})=\mathbf{g}_{d_{1}}^{\top}\mathrm{Diag}(\mathbf{p})^{- 1}\mathbf{g}_{d_{2}}=\sum_{t}TF(t,d_{1})TF(t,d_{2})p_{t}^{-1}\,.$ $$

Let us also calculate the IDF for this models. Let us assume that terms are rare enough ( $p_{t}|d|\ll 1$). Then we have that,

$$ $\displaystyle DF(d)$ $\displaystyle=1-\Pr(t\notin d)=1-(1-p_{t})^{|d|}\approx|d|p_{t}\,.$ $$

and therefore up to a scaling factor, the above expression is approximately equal to the $TF\cdot IDF$ relevance in (14).

## Appendix E More examples for comparison of LiSSA and PBRF

Here we present a complete list of 25 train images for comparison of PBRF and LiSSA. We calculate the LiSSA according to the hyperparameter recommendations in Table 1. Figure 5 shows scatter plots of LiSSA and PBRF for ResNet-18, and Figure 6 shows scatter plots for ResNet-50. Figure 7 shows reference images from the ImageNet dataset. Note that for one train image in Figure 6 the LiSSA got float overflow, we attribute it to the high value of the gradient norm.

Figure: Figure 5: Comparison of PBRF and LiSSA influence on ResNet-18 for 25 random train images. Each graph shows influence of one train image w.r.t. to 500 other test images. Reference number is show above the image, refer to Figure 7. The results are for ResNet-18, the $x$-axis is LiSSA, and the $y$-axis is PBRF.
Refer to caption: extracted/5879605/figures/retrain_25_resnet18.png

Figure: Figure 6: Same as Figure 5, but for ResNet-50.
Refer to caption: extracted/5879605/figures/retrain_25_resnet50.png

Figure: Figure 7: Reference images from ImageNet. Numbers above each image corresponde to those Figures 5, 6.
Refer to caption: extracted/5879605/figures/retrain_examples.png

## Peer Reviews

confidence: 3
contribution: 2
initial_rating: 5
presentation: 3
questions: - In section 5, the paper claims that the idea of inverse Hessian makes sense because small eigenvalues correspond to more specific and informative content. However, in traditional interpretation, small eigenvalues correspond to noise and carry little to no information. What is it not the case here? Moreover, are there references for the claim that large eigenvalues correspond to general language structure?
- The Hutchinson estimation and randomized sketching are used to estimate Tr(H) and lambda_max. In the experiments, how many Hessian-vector products are needed to get a good estimation?
- In Figure 3, lambda=5 is used, which seems to be large given that the GNH matrix in (7) is already SPSD. How does the results change using a smaller lambda?
soundness: 2
strengths: - The paper is well-organized and clearly written
- This paper gives the conditions on the hyperparameters for the LiSSA algorithm to converge. Such results are missing in the current literature.
summary: Influence function captures how the optimal model weights change with respect to a newly added data point and has the form of an inverse Hessian-vector product. LiSSA is classical solver to compute the influence function which essentially applying stochastic gradient descent to a quadratic objective. This paper proposes conditions to choose the hyperparameters (proximal hyperparameter, batch size, step size) of the LiSSA algorithm for computing the influence functions of a model to a newly added data point. The conditions guarantee the convergence of LiSSA in expectation. The convergence analysis shows that the batch size has to be sufficiently large in order for the LiSSA algorithm to converge.
weaknesses: - The contribution is somewhat limited because the convergence analysis seems to be standard, and the algorithm is only for computing the influence function.
- Some results in section 4 (Figure 1) seem to indicate that, in some situations, the influence function is not a good indicator of how the sensitive the optimal weights are with respect to newly added data points, as the LiSSA influence and PBRF influences are not correlated at all.
- The paper can be strengthened by reviewing more related work. E.g., are there other related algorithms for computing the influence function? Are there other variants of the influence function? Other work that analyze inverse Hessian-vector products (which are listed in the conclusion)? This can show the paper's position within the broader literature

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 2
questions: **Issues with the theoretical results**
- The submission does not address the large optimization literature on stochastic gradient descent, which has derived many similar convergence guarantees. The rates are not typically written in terms of a linear + constant term, as this is not convergent, and the literature is typically concerned with decreasing step-sizes (see e.g. [Bach and Moulines, 2013](https://arxiv.org/pdf/1306.2119)) or other schemes that can guarantee convergence. But the linear + constant result still appears regularly, for example see Theorem 2.1 in the work of [Needell et al., 2015](https://arxiv.org/pdf/1310.5715), Theorem 1 in the work of [Dieuleveut et al., 2016](https://arxiv.org/pdf/1602.05419), of Theorem 5.8 in the handbook of [Garrigos and Gower](https://arxiv.org/pdf/2301.11235). Given those results, the additional insights provided by the current submission appears small, but I might have missed a point that is specific to the inverse Hessian-product setting. If so, a discussion of existing literature in optimization and highlighting what makes the iHVP problem special would make the submission stronger.

- The text uses the word "converges" loosely, e.g. in corollary 1, "the algorithm converges in $T = \Omega(1/(\eta\lambda))$ steps". Taken literally, this is wrong, as the algorithm is not guaranteed to return $u^t = u^*$ after $T$ steps. I assume this is meant to say that the approximation error is bounded by some quantity after $T$ steps, but it isn't clear from the current writing. Changing the statement to a typical statement of "to achieve an error of $\frac{1}{2}\| u^t - u^*\|^2 \leq ...$, the algorithm should be run with parameters $\eta, B = ...$ for $T \geq ...$ iterations" would make it clearer. If the number of iterations is an upper bound, tt should also be $O$ rather than $\Omega$.

- The condition C.1 seems convenient but it isn't clear from the main text why it is introduced, or why it should hold. Extracting the key argument from the appendix to the main text to provide an explanation for the condition and why it is expected to hold would help readers understand the contribution. That the covariance scales with the trace appears to imply/assume that the covariance and the Hessian are the same, as in the toy quadratic model of [Zhang et al.](https://proceedings.neurips.cc/paper/2019/file/e0eacd983971634327ae1819ea8b6214-Paper.pdf) (see §3.5), and a citation to other people using similar assumptions might help make the case.

- The text "confirms the inverse scaling with batch size". Since the LHS of (C.1) is a variance term, the scaling as $\propto 1/B$ does not seem to need confirmation. The assumption that it is scaling with the trace however does seem to warrant empirical evidence. The scaling with the trace of the Hessian is the less intuitive part of the condition, and adding evidence for the scaling would strengthen the claim. 

**Writing**
The writing would benefit from an additional pass to improve the clarity of the message. Some of the statements appear wrong or misleading and the text contain typos or broken sentences. I only list a few examples below.

- Re: "Koh and Liang (2017) proposed a stochastic iterative approach called [...] (LiSSA, (Agarawal et al. 2017))". A more accurate description would be that Koh and Liang propose to _use_ the LiSSA algorithm of Agarawal et al., rather than propose the algorithm.
- Description such as "Basu et al. (2020) criticize the LiSSA [algorithm], suggesting that it lacks convergence for deep networks" and "we carefully analyze the convergence of the Lissa [algorithm]" should be made more specific. It is not clear whether this is refering to theoretical properties (will the algorithm actually converge to the solution if run sufficiently long) or whether the term "convergence" is used to refer to empirical performance, closer to "suggesting that the method performs poorly on deep networks" and "we analyze the [empirical performance] of the LiSSA algorithm as a function of its hyperparameters". 
- The sentence "SGD is known to work at least as well as the full gradient descent, even in terms of the number of updates (Harvey et al., 2019)" is misleading. Harvey et al. treat the non-smooth convex case, which is significantly different from the quadratic case considered here. The sentence should be removed.
- The equality in Eq. (3) should be replaced by an $\approx$ as it is the result from a truncated Taylor expansion
- L27-28, "interpret the internal computations in an understandable to a human way"
- "Generally speaking this is different from the empirical FIM [...], however, for low noise distributions the two might be used interchangeably (Martens, 2020)". I do not recall Martens making this argument, and it is not clear what is meant by "low-noise distribution".
soundness: 2
strengths: That the convergence theory for SGD can automate hyperparameter selection is a nice observation that is likely to be helpful to the community.
summary: The submission presents an analysis of the convergence of SGD for solving the inverse Hessian-vector-product, useful in the context of computing influence functions. The submission presents convergence results that depend on statistics of the Hessian and the batch size, and show that using those to set hyperparameters to achieve a desired accuracy leads to good performance, automating the choice of hyperparameters.
weaknesses: The main weakness is the first contribution of the submission, the convergence analysis in section 3. The presentation of the results ignores previous work on stochastic optimization and the description of the results lacks clarity. The main novelty would be the boxed condition C.1. However, it is unclear why this bound holds, as there is no theoretical explanation nor empirical validation provided in the main text. See detailed comments in Questions below.

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: Since this topic is about hyerparameter tuning, can we connect the proposed theory/results with iterative gradient-based optimization methods? More specifically, I know that for optimization algorithms like GD, its behavior will be influenced by the Hessian near equilibrium (i.e., the condition number of Hessian). So I am curious if this phenomenon is related to what is proposed in this paper.
soundness: 3
strengths: - The authors show that the two specific spectral properties of the Hessian, namely its trace and largest eigenvalue, determines the hyper parameters like the scaling factor, the batch size and the number of steps. So they evaluate these two quantities using random sketching. The conclusion is that batch size has to be sufficiently large for the LiSSA to converge. 

- code is open-source shared
summary: This paper studies influence functions, which are a popular tool for attributing models' outputs to training data. This reduces to calculation of inverse Hessian-vector products (iHVP) using the classical "Linear Time Stochastic Second-order Algorithm" (LiSSA) method. The authors show that the two specific spectral properties of the Hessian, namely its trace and largest eigenvalue, determines the hyper parameters like the scaling factor, the batch size and the number of steps. So they evaluate these two quantities using random sketching. The conclusion is that batch size has to be sufficiently large for the LiSSA to converge.
weaknesses: Overall I think the paper looks good.

---

confidence: 4
contribution: 2
initial_rating: 6
presentation: 3
questions: Regarding the weaknesses:
Weakness 1): There's a lot of batch-SGD analyses which at first glance may already contain the theorem.  For example:
A) Jain, Kakade, Kidambi, Netrapalli, Sidford. "Parallelizing.." 2018 JMLR
B) Varre, Pillaud-Vivien, Flammarion. "Last iterate convergnece.." 2021 NeurIPS
C) Ma, Bassily, Belkin.  "The power of Interpolation..." 2018 ICML
Can you argue why your analysis is novel, comparing to existing literature like the above?

Weakness 2): Can you show why the theorems apply to the setups in B3? Currently, this is `almost' done but not completely.  In particular the sentence on (l243): ``we find this condition to hold for classification with independent sampling" appears unfulfilled."

Weakness 3): Can you show empirically that C1 occurs, in the norm sense?  (tracial vs norm-type behavior is rarely equivalent in random matrix contexts).

Weakness 4a): Figure 1 is appears to be discussed in the text (l368-370) and 3 cases are explained. Considering Figure 1, why are these the 3 cases, (especially, which illustrates case '3'? Where do you argue that 'PBO _finetunning_ stirs away the model too far for the quadratic approximation to hold'. (and what does that mean?)

Weakness 4b): Figure 3.  Why does the difference picture show what you claim it does?  (In particular why do you average over the 10x10 square, and why are you not measuring the diagonals?). When you write (l464-466). ``As a result the influence similarity between an orignal sentence and a rewritten one appears to be consistently higher than between unrelated sentence."  Can you explain this?

Other questions: 
1) Can you explain precisely how you picked the 'batch' in Table 1, and moreover, what is the general procedure you suggest for selecting LiSSA parameters?  (It seems you _almost_ do this, but being extra-clear would help, especially as regard the batch.)

Minor:
1) (l80-83). ``Influence functions are calculated under the assumption that the optimized parameter $\theta$ of the model delivers minimum to the training loss."  This sentence is mangled. 
2) l409-410.  ``the choice of _dumping_ parameter''
3) l424: ``For example, a plausible interpretation of the directions $v_j$ would be that the top directions correspond to general language coherence and... small eigenvalues could correspond to more specific, informative content."  While the statement about large eigenvalues makes sense, the complementary statement about small eigenvalues does not -- first most small eigenvalues are not interpretable, and are the result of either diffuse (non-one-dimensional) embeddings or undesirable feature distortions.  Perhaps it would be better to just say you can screen out the large eigenvalues which are more about general language coherence.
soundness: 3
strengths: The paper does well to describe the end-to-end process of using inverse-Hessian influence metric on LLMs and Resnets, for which there are substantial empirical and theoretical questions to address.  It would appear that the cases that the authors considered are evaluated in reasonable level of detail, and that the empirical aspects of the paper should be reproducible.

While the main contributions are about methodology of using inverse-Hessian influences; section 5 provides some additional scientific value on the potential value of using inverse-Hessian influences.  While it is more speculative in nature, it provides some problem setups which might be the setting for future work.

The paper is also well-integrated into existing theory on influence functions, which brings additional value to the paper.

The main theorem seems to provide a suggested starting point for how to pick the parameters in LISSA (with explicit parameters).
summary: The paper considers inverse Hessian-vector products for the purpose of evaluating the influence of training data $(x_m, y_m)$, a data-label pair, on fitted parameters $\theta^*. 

For application to empirical problems, in particular language models, the Hessian is replaced by the Gauss-Newton matrix, and a regularized inverse-Hessian product is considered $(\lambda + H)^{-1}$. Extensive details are provided on the numerical considerations taken in evaluating Hessian-vector products.

An algorithm, LISSA, for estimating the Hessian-vector-product on a training point is considered.  It is effectively a version of minibatch SGD on a quadratic problem (8), where the gradient estimator is specific to this problem instance.  A convergence analysis in expectation and norm-convergence (up to a neighborhood) is provided.  These are taken under a relatively strong norm-type assumption on the minibatch gradients (C.1). 

This norm-type assumption is discussed in theory in the appendix in two settings and also empirically.

A table of hyperparameters are provided for a range of reasonable setups (Table 1), which (appear to be) conforming to the assumption C.1.  These require estimating the trace and operator norm of the hessian, and some discussion of the sketching procedure taken to do this is provided.

The convergence of LISSA is tested with these hyperparameters, to show good performance, and the influence is measured against a second method (PBRF).

Finally some discussion of the value of inverse-Hessian influence (vs gradient-based influence measures) are provided in Section 5.
weaknesses: I find the following overall weaknesses: 

1) The main theoretical contribution of the paper is a convergence analysis for the mini-batch SGD derived.  This is done under a very strong hypothesis, which is effectively a contractivity assumption in norm (C.1) or (9) in Theorem 1.  This makes the main theorem relatively easy to prove, but it also begs the question of showing when this satisfied.   Separately, mini-batch SGD on quadratics are studied to death, and it will be quite difficult to argue that this is truly a new analysis

2) Regarding the contractivity assumption: it is shown that in some sense the assumptions of Theorem 1 necessary, but this sense is effectively adversarial in the test-direction u (B.4Lemma 1).     On the other hand, the positive cases (B3 Motivations for C.1) do not show that C1 is satisfied.  It would substantially improve the paper if the theory covered these two cases.  (If it does, this should be explained).

3) The empirical results of the contractivity (Figure 4) do not consider C.1 but rather the trace.  Trace type conditions are much easier to satisfy, and unfortunately give little insight whether C.1 in fact holds.  A much more desirable theorem (to my taste -- and also more consistent with the experiments) would assume trace-type conditions on H and H^2, and some randomness assumptions on u that lead to the claimed convergence.  There is a lingering sense that C.1 and (9) of Theorem 1 is unverifiable.

4) The empirical conclusions in Section 4 and 5 are a little bit murky, and it's hard to align what appears in the simulations with what appear in the text.

Meta-review: This paper analyzes the convergence for LiSSA algorithm (used to compute inverse Hessian-vector products, which are useful for computing influence functions). The results show that a large (but mild) batch is required for convergence, and the optimal hyperparameters for the algorithm depends on spectral properties of the Hessian  (trace and largest eigenvalue). The paper then proposes to estimate these values by sketching method and give a practical way of finding the hyperparameters. The reviewers generally agree that this way of finding hyperparameters is interesting and practical. There are some concerns about whether this is only for a specific problem of estimating the influence function, or the novelty of the theoretical results. The author response clarified many aspects of the paper. However a couple reviewers are still not convinced.

Decision: Reject