# Generalized Leverage Score for Scalable Assessment of Privacy Vulnerability

Valentin Dorseuil, Jamal Atif, Olivier Cappé

► To cite this version:

Valentin Dorseuil, Jamal Atif, Olivier Cappé. Generalized Leverage Score for Scalable Assessment of Privacy Vulnerability. 2026. hal-05513289

HAL Id: hal-05513289

https://hal.science/hal-05513289v1

Preprint submitted on 16 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Generalized Leverage Score for Scalable Assessment of Privacy Vulnerability

Valentin Dorseuil*¹, Jamal Atif², and Olivier Cappé¹

¹DI ENS, École normale supérieure, Université PSL, CNRS, 75005 Paris, France

²CMAP, École polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France

## Abstract

Can the privacy vulnerability of individual data points be assessed without retraining models or explicitly simulating attacks? We answer affirmatively by showing that exposure to membership inference attack (MIA) is fundamentally governed by a data point's influence on the learned model. We formalize this in the linear setting by establishing a theoretical correspondence between individual MIA risk and the leverage score, identifying it as a principled metric for vulnerability. This characterization explains how data-dependent sensitivity translates into exposure, without the computational burden of training shadow models. Building on this, we propose a computationally efficient generalization of the leverage score for deep learning. Empirical evaluations confirm a strong correlation between the proposed score and MIA success, validating this metric as a practical surrogate for individual privacy risk assessment.

## 1 Introduction

Modern machine learning models, and deep neural networks in particular, are known to memorize aspects of their training data (Zhang et al., 2017; Carlini et al., 2019). This memorization induces privacy vulnerabilities that can be exploited by Membership Inference Attacks (MIAs), which aim to determine whether a specific data point was included in the training set (Shokri et al., 2017; Carlini et al., 2022). A principled defense against membership inference is provided by Differential Privacy (DP) (Dwork, 2006), implemented in deep learning through noise-injected stochastic gradient methods (Abadi et al., 2016). However, controlling the trade-off between privacy protection and model utility remains challenging. Noise calibration

typically relies on worst-case theoretical accounting, paired with empirical privacy auditing via MIAs, and often leads to either over-conservative noise levels or insufficient privacy protection. In this context, designing MIAs for empirical auditing is essential to quantify leakage in non-private models (Yeom et al., 2018) or to validate the practical tightness of DP guarantees in private models (Nasr et al., 2021; Jagielski et al., 2020).

While such auditing is now standard practice (Carlini et al., 2022; Nasr et al., 2021; Zarifzadeh et al., 2024), relying on aggregate metrics like average accuracy or AUC is insufficient. Such global measures can obscure critical risk heterogeneity, as outliers and rare subgroups are significantly more prone to memorization than typical samples (Carlini et al., 2022; Feldman and Zhang, 2020). Consequently, a model certified as private on average may still expose specific points to high privacy risks.

To address this heterogeneity, recent work is increasingly focusing on *individual* privacy risk assessment, aiming to quantify membership leakage at the level of each data point, separately, rather than in aggregate. State-of-the-art methods for per-sample auditing (Carlini et al., 2022; Zarifzadeh et al., 2024) predominantly rely on *shadow models*. These techniques train multiple reference models on random data splits to characterize each data point's influence on the model's behavior. This enables identification of which points exhibit increased sensitivity to their presence in the training dataset. However, these techniques are computationally prohibitive, especially for large-scale models, as they require retraining the model multiple times. This leads to our central question stated in the abstract: *can the privacy vulnerability of individual data points be assessed without retraining models or explicitly simulating attacks?*

In classical statistics, the *leverage score* quantifies a data point's geometric influence on a model, independent of its label. We establish that in Gaussian linear

<footnote>*Corresponding author: valentin.dorseuil@ens.psl.edu
Preprint. Under review.</footnote>
models, this score precisely characterizes membership
inference vulnerability. We prove that under an op-
timal black-box attack, the privacy loss distribution
is controlled by a single scalar, the leverage score. In
other words, membership inference vulnerability is
fundamentally about self-influence: samples that are
geometrically positioned to have a disproportionate
effect on the model's learned parameters are the most
at risk of privacy leakage. Consequently, while spe-
cific outcomes fluctuate due to noise, the individual
average privacy risk in the linear regime is determined
by data geometry.

To extend this analysis to deep neural networks,
we introduce the *Generalized Leverage Score* (GLS).
Derived via implicit differentiation of the training opti-
mality conditions, the GLS measures the infinitesimal
sensitivity of a model’s prediction to its own label,
generalizing the leverage score to both regression and
classification settings. While the exact computation
is expensive for deep networks, we demonstrate that
a last-layer approximation remains highly effective
in practice. This allows us to compute a scalable,
theoretically principled, proxy for privacy risk that
correlates with the success of state-of-the-art attacks,
without the need for retraining or shadow models.

**Contributions.** (i) We prove that for Gaussian linear models under black-box access, the leverage score is the sufficient statistic characterizing both the privacy loss distribution and the optimal membership inference test. (ii) We extend this leverage score to general differentiable models via the *Generalized Leverage Score* (GLS), deriving a principled and scalable estimator of privacy vulnerability. (iii) Through multiple experiments, we show that this metric serves as a surrogate for individual privacy risk. It identifies vulnerable samples, showing a strong correlation with state-of-the-art shadow model attacks, at a reduced computational cost.

2 Related Work

**Membership Inference Attacks (MIA).** Membership inference aims to determine if a specific sample was used to train a model. Early approaches relied on simple metric-based classifiers, exploiting overfitting signals such as prediction confidence, entropy, or the magnitude of gradients (Shokri et al., 2017; Yeom et al., 2018). While computationally inexpensive, these methods often struggle to distinguish between a “vulnerable” member and a “hard” non-member.

To address this, current state-of-the-art methods
adopt the *shadow model* paradigm Shokri et al. (2017).

By training multiple models on different data splits,
attacks like LiRA (Likelihood Ratio Attack) (Carlini
et al., 2022) essentially perform a hypothesis test on
the loss distribution of a target point (Zarifzadeh
et al., 2024). While highly effective, this approach is
computationally prohibitive, often requiring hundreds
of training runs to estimate the risk for a single point
accurately. Our work seeks to achieve the precision
of these hypothesis-test based methods without their
computational burden, by substituting retraining with
geometric analysis.

Influence Functions and Data Attribution.

Originating in robust statistics with influence mea-
sures such as Cook's distance (Cook, 1977), influence
functions quantify the effect of a training point on
model parameters or predictions. In modern deep
learning, Koh and Liang (2017) reintroduced influ-
ence functions via Hessian-vector products to explain
model's behavior and identify mislabeled data.

Recent work has begun to explore the connection
between influence and privacy. Notably, Feldman and
Zhang (2020) and Feldman (2020) link the memoriza-
tion of a sample to its influence on the learning process,
arguing that memorization is necessary for general-
ization in long-tailed distributions. However, classical
influence functions typically measure the effect of a
training point on a separate test point. In contrast,
our GLS formulation focuses on self-influence (the
sensitivity of a prediction to its own label), identi-
fying this specific form of leverage as the canonical
driver of privacy leakage.

**Privacy Auditing and Heterogeneity.** Differen-
tial Privacy (DP) (Dwork, 2006) provides worst-case
guarantees for membership privacy. Using algorithms
such as DP-SGD (Abadi et al., 2016), one can train
Deep Learning models to be private. However, these
bounds are often loose and do not capture the empir-
ical reality that privacy risk is non-uniform: some
outlier samples are far more exposed than others
(Zarifzadeh et al., 2024).

Auditing individual privacy risk (estimating the
specific (ε, δ) for a given sample) remains an open
challenge. Existing auditing tools rely heavily on
the aforementioned shadow model techniques or ran-
domized smoothing (Lecuyer et al., 2019), which are
difficult to scale. Our work contributes to this domain
by proposing an individual privacy risk metric for
identifying vulnerable samples efficiently. We believe
that this approach can complement existing audit-
ing methods, providing a first-pass filter to identify
high-risk points without retraining.
# 3 Leverage Scores and Membership Vulnerability

We begin by analyzing how data heterogeneity impacts the effectiveness of Membership Inference Attacks (MIAs) in the fixed-design Gaussian linear regression framework. This setting is selected for its analytical tractability and its ability to represent data variance through the fixed design matrix. In contrast to approaches that assume an i.i.d. data distribution, our fixed-design analysis captures the non-uniform risk exposure inherent to heterogeneous datasets.

We provide a characterization of the optimal membership inference test and its trade-off curve in this setting, where the average is taken over all possible training noise realizations for a fixed dataset, rather than across data points. By quantifying vulnerability using leverage scores, we establish a rigorous link between a point’s structural influence and its individual privacy risk exposure.

## 3.1 Problem Setting and Definitions

Consider a fixed design matrix $X \in \mathbb{R}^{n \times d}$, where each row $x_i^\top$ represents a data point. We examine the multivariate linear model:

$$Y = X\Theta^* + E, \quad (1)$$

where $Y \in \mathbb{R}^{n \times m}$ is the response matrix, $\Theta^* \in \mathbb{R}^{d \times m}$ is the true parameter matrix, and $E \in \mathbb{R}^{n \times m}$ represents noise. We assume the noise is centered, independent and identically distributed (i.i.d.) Gaussian, such that for any row $i$, $E[E_{i.}] = 0$ and $Cov(E_{i.}) = \sigma^2 I_m$. We assume $n \ge d$ and that $X$ has full column rank.

The Ordinary Least Squares (OLS) estimator is $\hat{\Theta} = (X^\top X)^{-1} X^\top Y$. The fitted values are $\hat{Y} = HY$, where $H = X(X^\top X)^{-1}X^\top$ is the hat matrix. For a specific data point $i$, the residual vector $r_i \in \mathbb{R}^m$ is defined as:

$$r_i = y_i - \hat{\Theta}^\top x_i = \sum_{j=1}^{n} (\delta_{ij} - h_{ij}) y_j, \quad (2)$$

where $\delta_{ij}$ is the Kronecker delta. The diagonal elements of the hat matrix, $h_{ii} = x_i^\top(X^\top X)^{-1}x_i$, are the leverage scores. They satisfy $0 \le h_{ii} \le 1$ and $\sum_{i=1}^n h_{ii} = d$, serving as a measure of the geometric influence of $x_i$ on the model’s predictions (Belsley et al., 1980).

## 3.2 Distinguishing Members from Non-Members

To formalize the Membership Inference Attack, we define two hypotheses for a data point $(x_i, y_i)$. Un-

der the member hypothesis ($H_1 = H_{train}$), the pair was included in the training set used to compute $\hat{\Theta}$. Under the non-member hypothesis ($H_0 = H_{test}$), the model was trained on the dataset, but we evaluate the residual on an independent observation $\tilde{y}_i$ generated from the same distribution: $\tilde{y}_i = \Theta^{*\top}x_i + \varepsilon'_i$, where $\varepsilon'_i$ is independent noise.

**Proposition 3.1 (Residual Distributions).** Under the Gaussian noise assumption, the residuals for the i-th data point under the member and non-member hypotheses follow distinct multivariate normal distributions:

$$r_i | H_{train} \sim N(0, \sigma^2(1 - h_{ii})I_m), \quad (3)$$

$$r_i | H_{test} \sim N(0, \sigma^2(1 + h_{ii})I_m). \quad (4)$$

Consequently, the squared norms of the residuals follow scaled Chi-squared distributions:

$$\|r_i\|^2 | H_{train} \sim \sigma^2(1 - h_{ii})\chi^2(m), \quad (5)$$

$$\|r_i\|^2 | H_{test} \sim \sigma^2(1 + h_{ii})\chi^2(m). \quad (6)$$

The proof is provided in Appendix A.2. This result highlights that members have lower residual variance than non-members, with the gap controlled explicitly by the leverage score $h_{ii}$. In particular, the distributions differ in scale and not in location, leading to fundamentally non-symmetric error trade-offs in membership inference. This asymmetry leads to easier detection of non-members (low false negative rate $\beta$) than members (low false positive rate $\alpha$), which is consistent with the empirical behavior of state-of-the art MIAs (Carlini et al., 2022).

## 3.3 A Membership Inference Attack Perspective

Using Proposition 3.1, we derive the theoretically optimal membership inference attack. The Neyman-Pearson Lemma (Neyman and Pearson, 1933) states that the most powerful test for distinguishing these hypotheses is the Likelihood Ratio Test.

Due to the heterogeneity of the data (varying $h_{ii}$), a global threshold on the raw loss $\|r_i\|^2$ is suboptimal. Instead, the optimal decision rule requires a sample-specific normalization.

**Proposition 3.2 (Optimal MIA Test).** Let $S_i(r_i)$ be the log-likelihood ratio statistic for the i-th sample. The most powerful test at level $\alpha$ rejects the non-member hypothesis $H_{test}$ (predicts membership) if $S_i(r_i) > \gamma_\alpha$, where the sufficient statistic is given by:

$$S_i(r_i) = \frac{m}{2} \ln \left( \frac{1 + h_{ii}}{1 - h_{ii}} \right) - \frac{h_{ii}}{\sigma^2(1 - h_{ii}^2)} \|r_i\|^2. \quad (7)$$
The proof of this result is provided in Appendix A.3. To obtain the optimal $(\alpha, \beta)$ trade-off curve, samples must be ranked in descending order of their score $S_i(\mathbf{r}_i)$. This score represents a data-dependent affine transformation of the loss $\|\mathbf{r}_i\|^2$. This relationship aligns with the design of parametric MIAs that fit Gaussian distributions to scores. It suggests that such methods implicitly learn this affine transformation to distinguish members from non-members (Carlini et al., 2022; Zarifzadeh et al., 2024).

We follow the approach of f-Differential Privacy defined in Dong et al. (2022) to characterize the membership inference trade-off curve as a functional relationship between error types.

**Proposition 3.3 (MIA Errors Curve).** For the point $i$, the trade-off between the false alarm rate $\alpha$ (Type-I error) and the missed detection rate $\beta$ (Type-II error) is described by the curve:

$$ \beta_i(\alpha_i) = 1 - F_m \left( \frac{1 + h_{ii}}{1 - h_{ii}} F_m^{-1}(\alpha_i) \right). \quad (8) $$

where $F_m$ is the Cumulative Distribution Function of the $\chi^2(m)$ distribution.

This result highlights that $h_{ii}$ is the unique parameter controlling the distinguishability of training points from test points in the Black-Box setting, i.e., when the attacker only has access to the loss values. The proof is provided in Appendix A.4 with an illustration for multiple $h_{ii}$ values in Figure 5.

## 3.4 An Influence Point of View

The vulnerability of high-leverage points established in Proposition 3.2 can be intuitively understood through the lens of self-influence.

While the leverage score $h_{ii}$ was defined geometrically via the projection matrix, it admits an operational interpretation as the self-influence of the $i$-th observation.

With the OLS estimator, $\hat{\Theta}$ we have:

$$ \hat{\mathbf{y}}_i = \sum_{j=1}^{n} h_{ij} \mathbf{y}_j = \sum_{j=1}^{n} \mathbf{x}_i^{\top} (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{x}_j \mathbf{y}_j. \quad (9) $$

The fitted value for point $i$ is a linear combination of all observed labels $\mathbf{y}_j$. Thus leverage score $h_{ii}$ corresponds to the sensitivity of the $i$-th fitted value $\hat{\mathbf{y}}_i$ with respect to the observed target $\mathbf{y}_i$. This yields the identity:

$$ \frac{\partial \hat{\mathbf{y}}_i}{\partial \mathbf{y}_i} = h_{ii} \mathbf{I}_m. \quad (10) $$

It quantifies how much the model's prediction for a specific point changes when the label of that point

is perturbed. The scalar leverage score can then be recovered via the trace: $h_{ii} = \frac{1}{m} \text{Tr} \left( \frac{\partial \hat{\mathbf{y}}_i}{\partial \mathbf{y}_i} \right)$.

This identity highlights the link between the leverage score and the model's sensitivity to its own training labels. For high-leverage points ($h_{ii} \approx 1$), the model is forced to interpolate the observation $\mathbf{y}_i$ to minimize the global squared error. The specific noise instance $\varepsilon_i$ present in a member's label is directly encoded into the prediction $\hat{\mathbf{y}}_i$. Consequently the training residual $\mathbf{r}_i$ vanishes, making membership easily detectable.

While the geometric definition of $h_{ii}$ is specifically relevant to OLS, the sensitivity formulation in Equation (10) is general, allowing us to extend leverage scores to non-linear deep neural networks.

# 4 Generalized Leverage Scores (GLS)

We formalize the Generalized Leverage Score (GLS) to extend the classical sensitivity interpretation to arbitrary differentiable models.

## 4.1 Formal Framework and Definition

We consider a dataset $\mathcal{D} = \{(\mathbf{x}_i, \mathbf{y}_i)\}_{i=1}^n$ and a model $f_\theta: \mathcal{X} \to \mathcal{Y}$ parameterized by $\theta \in \Theta \subset \mathbb{R}^p$. The output space $\mathcal{Y}$ is $m$-dimensional (for example $m = 1$ for scalar regression and $m = C$ for $C$-class classification). Training is performed by minimizing the empirical risk:

$$ \mathcal{L}(\theta) = \frac{1}{n} \sum_{i=1}^{n} \ell(f_{\theta}(\mathbf{x}_i), \mathbf{y}_i), \quad (11) $$

where $\ell$ is a twice-differentiable loss function. Let $\hat{\theta}$ denote the optimal parameters obtained at convergence, i.e., $\hat{\theta} \in \arg \min_{\theta \in \Theta} \mathcal{L}(\theta)$.

**Definition 4.1 (Generalized Leverage Score).** The Generalized Leverage Score (GLS) of the $i$-th sample is defined as the infinitesimal sensitivity of the model's prediction at $\mathbf{x}_i$ to its own observed label $\mathbf{y}_i$:

$$ \text{GLS}_i = \frac{\partial f_{\theta}(\mathbf{x}_i)}{\partial \mathbf{y}_i} \in \mathbb{R}^{m \times m}. \quad (12) $$

This metric captures a self-influence effect; it quantifies how much the model would alter its prediction for a specific sample if that sample's label was perturbed, accounting for the implicit change in the learned parameters $\hat{\theta}$.
## 4.2 Analytic Derivation via Implicit Differentiation

Computing GLS$_i$ directly is challenging as $\hat{\theta}$ depends on $\mathbf{y}_i$ implicitly through the optimization process. To isolate this dependence, we follow the approach of Koh and Liang (2017) and derive a closed-form expression for GLS$_i$.

**Proposition 4.2 (Closed-form GLS).** Assuming the averaged loss function $\mathcal{L}$ is twice differentiable and its Hessian $\mathbf{H}_{\hat{\theta}}$ is invertible, the GLS for a univariate output is given by:

$$ \text{GLS}_i = -\frac{1}{n} \|\nabla_{\theta} f_{\hat{\theta}}(\mathbf{x}_i)\|_{\mathbf{H}_{\hat{\theta}}^{-1}}^2 \left(\frac{\partial^2 \ell}{\partial \mathbf{y} \partial f}\right), \quad (13) $$

where $\|\mathbf{v}\|_A^2 = \mathbf{v}^\top A \mathbf{v}$ denotes the squared Matrix norm. For multivariate outputs, this generalizes to:

$$ \text{GLS}_i = -\frac{1}{n} \mathbf{J}_{\theta} f_{\hat{\theta}}(\mathbf{x}_i) \mathbf{H}_{\hat{\theta}}^{-1} \mathbf{J}_{\theta} f_{\hat{\theta}}(\mathbf{x}_i)^{\top} \left( \frac{\partial^2 \ell}{\partial \mathbf{y} \partial f} \right), \quad (14) $$

where $\mathbf{J}_{\theta} f_{\hat{\theta}}(\mathbf{x}_i) \in \mathbb{R}^{m \times p}$ is the Jacobian of the model outputs with respect to the parameters, $\mathbf{H}_{\hat{\theta}}$ is the Hessian of the loss function $\mathcal{L}$, i.e. $\mathbf{H}_{\hat{\theta}} = \nabla_{\theta}^2 \mathcal{L}(\hat{\theta})$, and $\ell$ denotes the pointwise loss evaluated at $(f_{\hat{\theta}}(\mathbf{x}_i), \mathbf{y}_i)$.

This result follows from applying implicit differentiation to the first-order optimality condition $\nabla_{\theta} \mathcal{L}(\hat{\theta}) = 0$, noting that $\hat{\theta}$ depends implicitly on $\mathbf{y}_i$ and $\mathcal{L}$ depends explicitly on $\mathbf{y}_i$ through $\ell(f_{\hat{\theta}}(\mathbf{x}_i), \mathbf{y}_i)$. The full derivation is provided in Appendix B.2.

The scaling factor $\frac{\partial^2 \ell}{\partial \mathbf{y} \partial f}$ acts as a loss-specific weighting term in the GLS expression. We derive this term for common loss functions in both regression and classification settings.

## 4.3 GLS for Regression Models

In the multivariate regression setting, we typically consider the quadratic loss. For a single data point $(\mathbf{x}, \mathbf{y})$ the loss is defined as: $\ell(f(\mathbf{x}), \mathbf{y}) = \|f(\mathbf{x}) - \mathbf{y}\|_2^2$.

**Proposition 4.3 (Multivariate Regression).** For the quadratic loss with targets $\mathbf{y} \in \mathbb{R}^m$ and predictions $f \in \mathbb{R}^m$, the second derivative is:

$$ \frac{\partial^2 \ell}{\partial \mathbf{y} \partial f} = -2\mathbf{I}_m. \quad (15) $$

This gives us:

$$ \text{GLS}_i = \frac{2}{n} \mathbf{J}_{\theta} f_{\hat{\theta}}(\mathbf{x}_i) \mathbf{H}_{\hat{\theta}}^{-1} \mathbf{J}_{\theta} f_{\hat{\theta}}(\mathbf{x}_i)^{\top}. \quad (16) $$

## 4.4 GLS for Classification Models

For classification, we consider the cross-entropy loss $\ell(\hat{\mathbf{p}}, \mathbf{y}) = -\sum_{j=1}^m y_j \log(\hat{p}_j)$, where $\mathbf{y} \in \{0, 1\}^m$ is the one-hot label vector and $\hat{\mathbf{p}} = s(f_{\hat{\theta}}(\mathbf{x})) \in [0, 1]^m$ is the vector of predicted probabilities obtained via the softmax function $s$.

The GLS depends on the loss through the second derivative $\frac{\partial^2 \ell}{\partial \mathbf{y} \partial f}$, which for cross-entropy reduces to a constant matrix:

**Proposition 4.4 (Cross-Entropy Derivatives).** For the cross-entropy loss with one-hot labels $\mathbf{y}$ and logits $f$, the second derivative is:

$$ \frac{\partial^2 \ell}{\partial \mathbf{y} \partial f} = -\mathbf{I}_m. \quad (17) $$

While the GLS is formally derived in terms of logit sensitivity, interpreting these scores in the probability space is often more practical. For classification, we thus consider the GLS in the probability space. Applying the chain rule, we have:

$$ \begin{aligned} \text{GLS}_i &= \frac{\partial \hat{\mathbf{p}}_i}{\partial \mathbf{y}_i} \\ &= \frac{\partial s}{\partial f} \times \frac{\partial f_{\hat{\theta}}(\mathbf{x}_i)}{\partial \mathbf{y}_i} \\ &= \frac{1}{n} \underbrace{(\text{diag}(\hat{\mathbf{p}}_i) - \hat{\mathbf{p}}_i \hat{\mathbf{p}}_i^\top)}_{\text{Softmax Jacobian}} \mathbf{J}_{\theta} f_{\hat{\theta}}(\mathbf{x}_i) \mathbf{H}_{\hat{\theta}}^{-1} \mathbf{J}_{\theta} f_{\hat{\theta}}(\mathbf{x}_i)^{\top}. \end{aligned} $$

This formulation yields a leverage score properly normalized within the probability simplex. We use this score as our primary privacy metric in the following experiments. Detailed derivations are in Appendix B.3.

## 4.5 Interpretation in Binary Classification

To further interpret the GLS in classification, we consider the classical binary logistic regression model. The model predicts $\hat{p}_i = \sigma(\hat{y}_i)$, where $\hat{y}_i = \mathbf{x}_i^\top \theta$ and $\sigma: z \mapsto 1/(1 + e^{-z})$ is the sigmoid function. The standard cross entropy loss is $\ell(\hat{y}_i, y_i) = -y_i \log(\hat{p}_i) - (1 - y_i) \log(1 - \hat{p}_i)$.

**Proposition 4.5 (Binary Logistic Leverage Score).** Under the cross entropy loss, for a binary logistic regression model, the GLS in the probability space reduces to:

$$ \frac{\partial \hat{p}_i}{\partial y_i} = \sigma'(\hat{y}_i) \times \text{GLS}_i = w_i \mathbf{x}_i^\top (\mathbf{X}^\top \mathbf{W} \mathbf{X})^{-1} \mathbf{x}_i, \quad (18) $$

with $w_i = \hat{p}_i(1 - \hat{p}_i) = \sigma'(\hat{y}_i)$ and $\mathbf{W} = \text{diag}(w_1, \dots, w_n)$.
The proof is provided in Appendix B.4. We recover here the leverage score from logistic regression defined in Pregibon (1981). This score can be interpreted as a weighted version of the input's norm in the Hessian-induced metric:

$$ \frac{\partial \hat{p}_i}{\partial y_i} = \sigma'(\hat{y}_i) \times \|\mathbf{x}_i\|_{H_{\hat{\theta}}^{-1}}^2, \quad (19) $$

with $H_{\hat{\theta}} = X^{\top}WX$ the Hessian of the loss at convergence.

The term $w_i = \hat{p}_i(1 - \hat{p}_i)$ reaches its maximum at $\hat{p}_i = 0.5$, which corresponds to sample located on the decision boundary. This aligns with the LiRA intuition that samples near the decision rule are most sensitive to perturbations in the training set, and thus most vulnerable to membership inference (Carlini et al., 2022).

# 5 Scalable Computation & Approximations

Computing the GLS directly from its closed-form expression can be computationally prohibitive for deep networks due to the need to invert the $p \times p$ Hessian matrix, where $p$ is the number of model parameters. To address this challenge, we propose an efficient algorithm that leverages Hessian-vector products (HVPs) and iterative solvers to compute the GLS without explicitly forming or inverting the Hessian.

## 5.1 Efficient GLS Computation

The key observation is that the GLS expression involves the term $H_{\hat{\theta}}^{-1}J_{\theta}f_{\hat{\theta}}(\mathbf{x}_i)^{\top}$. Instead of computing the inverse Hessian directly, we can reformulate this as solving a linear system:

$$ H_{\hat{\theta}} Z_i = J_{\theta} f_{\hat{\theta}}(\mathbf{x}_i)^{\top}, \quad (20) $$

for each sample $i$, where $Z_i$ is the unknown matrix we wish to compute. Once $Z_i$ is obtained, the GLS can be computed as:

$$ \text{GLS}_i = -\frac{1}{n} \left( \frac{\partial^2 \ell}{\partial \mathbf{y} \partial f} \right) J_{\theta} f_{\hat{\theta}}(\mathbf{x}_i) Z_i. \quad (21) $$

To solve the linear system efficiently, we employ the conjugate gradient (CG) method, which only requires the ability to compute Hessian-vector products. These products can be computed efficiently using automatic differentiation techniques without explicitly forming the Hessian matrix (Pearlmutter, 1994).

The overall algorithm for computing the GLS for a set of target samples $\mathcal{X} = \{\mathbf{x}_i\}_{i=1}^k$ is summarized

**Algorithm 1** Inverse-Hessian Jacobian Product via Conjugate Gradient Solver

**Input:** $\mathcal{D}_{train}$, Model $f_{\hat{\theta}}$, Jacobians $J^{\top}$, Damping $\lambda$, Iterations $T$

**Output:** $Z \in \mathbb{R}^{p \times m}$
Define the linear operator : $A[J^{\top}] = (H_{\hat{\theta}} + \lambda I)J^{\top}$
$Z \leftarrow \text{CG Algorithm}(\text{operator} = A, \text{target} = J^{\top}, \text{iters} = T)$

**return** $Z$

**Algorithm 2** Generalized Leverage Score Computation

**Input:** Target samples $\mathcal{X}$, $\mathcal{D}_{train}$, $f_{\hat{\theta}}$, Damping $\lambda$, Iterations $T$

**Output:** $\{\text{GLS}_i\}_{i=1}^k$
**for** $i = 1$ **to** $k$ **do**
&emsp;&emsp;$J_i \leftarrow \{\nabla_{\theta} f_{\hat{\theta}}(\mathbf{x}_i)\}_{i=1}^k$
&emsp;&emsp;$M_i \leftarrow -\frac{1}{n} \frac{\partial^2 \ell}{\partial \mathbf{y} \partial f}(f(\mathbf{x}_i), \mathbf{y}_i)$
&emsp;&emsp;$Z_i \leftarrow \text{Algorithm 1}(\mathcal{D}_{train}, f_{\hat{\theta}}, J_i^{\top}, \lambda, T)$
&emsp;&emsp;$\text{GLS}_i \leftarrow J_i Z_i M_i$

**end for**

**return** $\{\text{GLS}_i\}_{i=1}^k$

in Algorithm 2. The key steps involve computing the Jacobians for the target samples, solving the linear systems using CG (Algorithm 1), and finally assembling the GLS values. Using `torch.vmap` (Paszke et al., 2019) and tensor parallelization, we implemented a version of this algorithm that was able to process multiple target samples at the same time.

**Computational Complexity.** The proposed framework avoids the $O(p^3)$ time and $O(p^2)$ memory requirements of direct Hessian inversion by using a conjugate gradient (CG) solver. The primary computational cost is dominated by solving the linear system $A[J^{\top}]$, which computes Hessian-Matrix Products (HMPs) via a pass over the training set $\mathcal{D}_{train}$. For a model with $p$ parameters and $m$ output logs, the total time complexity is $O(kT \text{nm} \text{p})$, where $k$ is the number of target points, $n$ is the training set size, and $T$ is the number of CG iterations.

<footnote>The Hessian Jacobian product can be computed using mini-batches for the Hessian part, allowing for efficient estimation by summing the contributions from each batch without requiring the full dataset in one pass.</footnote>

<footnote>The CG algorithm iteratively refines the solution to the linear system using only matrix-vector products, making it suitable for large-scale problems where the matrix is not explicitly formed. The details of the CG algorithm can be found in Hestenes and Stiefel (1952).</footnote>
<table><caption>Table 1: Computational Complexity of GLS Subprocesses.</caption><thead><tr><th>Subprocess</th><th>Time</th><th>Memory</th></tr></thead><tbody><tr><td>Jacobian</td><td>O(pm)</td><td>O(pm)</td></tr><tr><td>Hessian Matrix Product</td><td>O(npm)</td><td>O(pm)</td></tr><tr><td>CG Solver (T iter)</td><td>O(Tnpm)</td><td>O(pm)</td></tr><tr><td><strong>Total for k points</strong></td><td>O(kTrpm)</td><td>O(pm)</td></tr></tbody></table>

## 5.2 Approximation via Layer Restriction

Exact computation of GLS over the full training set ($k = n$) is intractable for deep networks as our algorithm yields an $O(n^2)$ time complexity. Thus, our proposed algorithm is computationally feasible for a moderate number of target points $k$ since the time complexity scales linearly with $k$.

<figure><img src="image_2.png" /><figcaption>Figure 1: Spearman's correlation (green line) between LiRA scores and GLS (trace), and computational time (blue area) for different network depths. Error bars show 95% confidence intervals across 16 models. Computational times are measured on a single A100 GPU.</figcaption></figure>

To mitigate this, we explored computing our GLS by differentiating only through a subset of layers, and found that restricting the computation to the last few layers yields nearly the same correlation with our comparison metrics as the full model. Formally, this corresponds to computing $\frac{\partial y}{\partial y}$ on a restricted model $\theta_{sub}$ where the feature extractor is frozen. This is equivalent to masking the full Jacobian and Hessian to the subspace $\theta_{sub}$.

The results can be observed in Figure 1, where the correlation between the LiRA Score and the GLS is

plotted against computation time for various differentiation depths.

The resulting complexity is reduced to $O(kTnmp_{sub})$, where $p_{sub}$ is the number of parameters in $\theta_{sub}$.

In particular, when the last layer is linear, the GLS admits a closed-form solution. We recover the classical statistical leverage scores in the feature space, allowing for easier computation. With $d$ being the feature dimension, the closed-form lead to an improved complexity of $O(n d^2 + d^3)$ for the regression case with quadratic loss (even with $m > 1$) and $O(n m^2 d^2 + m^3 d^3)$ for the classification case with the cross entropy loss. More details can be found in Appendix C.

# 6 Experiments

In this section, we empirically evaluate the effectiveness of the Generalized Leverage Score (GLS) as a membership inference metric. We assess its correlation with the LiRA score, an established benchmark for membership inference attacks, and analyze its performance in identifying outlier samples within the training data.

## 6.1 Setup

**Dataset and Model.** We conduct experiments on the CIFAR-10 dataset (Krizhevsky and Hinton, 2009), using a ResNet-18 architecture (He et al., 2016) trained to achieve approximately 92% test accuracy. The training set consists of 50,000 images, while the test set contains 10,000 images.

In order to compute the LiRA scores, we follow the protocol established by Carlini et al. (2022), training an ensemble of 200 shadow models with architectures identical to the target model. Each shadow model is trained on a random subset of 40,000 images from CIFAR-10, with the remaining 10,000 images reserved for evaluation. The LiRA scores are then computed for all the samples, split between members and non-members of the target model's training set.

**Computation** To compute the Generalized Leverage Scores (GLS), we implement the algorithm outlined in Section 5.1. We set the damping factor $\lambda$ between $10^{-4}$ and $10^{-2}$, depending on the differentiation depth and run the conjugate gradient (CG) solver for $T = 100$ iterations. In most cases, the solver returned a solution (with a tolerance of $10^{-3}$ on the equality) before the maximum number of iterations

<footnote>The LiRA score is computed as the likelihood ratio of the observation under two Gaussian distributions modeling the member and non-member hypotheses. The parameters of these</footnote>

<footnote>distributions are estimated using an ensemble of shadow models. More information can be found in Carlini et al. (2022).</footnote>
<table><caption>Table 2: Spearman correlation between metrics computed using shadow models (LiRA) and Generalized Leverage Scores for different differentiation depths. Results are averaged over 16 different models; 95% quantiles are reported next to each entry using ± notation. Highest correlations per metric are in bold.</caption><thead><tr><th rowspan="2">Shadow Metrics</th><th colspan="4">Generalized Leverage Scores</th></tr><tr><th>last layer</th><th>2 last layers</th><th>4 last layers</th><th>full model</th></tr></thead><tbody><tr><td>LiRA Score</td><td>0.50 ± 0.05</td><td><strong>0.51</strong> ± 0.04</td><td><strong>0.51</strong> ± 0.03</td><td><strong>0.51</strong> ± 0.04</td></tr><tr><td>|σ<sub>test</sub>/σ<sub>train</sub>|</td><td>0.40 ± 0.04</td><td><strong>0.41</strong> ± 0.03</td><td><strong>0.41</strong> ± 0.04</td><td><strong>0.41</strong> ± 0.03</td></tr><tr><td>TPR@FPR = 0.05</td><td>0.27 ± 0.04</td><td>0.28 ± 0.04</td><td>0.28 ± 0.05</td><td><strong>0.29</strong> ± 0.05</td></tr><tr><td>|μ<sub>test</sub> - μ<sub>train</sub>|</td><td>0.06 ± 0.04</td><td>0.06 ± 0.03</td><td><strong>0.07</strong> ± 0.04</td><td>0.06 ± 0.04</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>Figure 2: t-SNE visualization of the model's final layer representations for the "Frog" class. Point coloring represents the leverage score magnitude, highlighting the geometric distribution of high-leverage samples.</figcaption></figure>

was reached. We trained 16 different models with different random seeds to compute the GLS scores and report averaged results. To compare results across models, we keep the first 128 samples in the train set for all models.

We tried different operations on the GLS matrix to get a scalar score: the Trace, the Frobenius Norm and the Spectral Norm. We observed that each of these operations lead to similar results. In the rest of the paper we report results with the Trace operator.

## 6.2 Metrics and Results

**Correlation with LiRA** Table 2 reports the Spearman correlation between the Generalized Leverage Scores, computed for different differentiation depth, and various shadow model indicators, including the LiRA score, the mean gap ($|\mu_{in} - \mu_{out}|$), the standard deviation ratio ($\sigma_{in}/\sigma_{out}$), and the TPR at low FPR. We observe a positive rank correlation between GLS and the LiRA score as well as other metrics. This confirms that the leverage score captures a similar vulnerability signal as more computationally expensive

shadow models: high-leverage points correspond to samples where attacks are more likely to confidently predict membership.

**Outlier Detection** To further understand the characteristics of high-leverage samples, we visualize the model's final layer representations using t-SNE (van der Maaten and Hinton, 2008). In Figure 2, we color the points based on their GLS magnitude, revealing that most of the low-leverage samples cluster around the class center in the representation space. In contrast the points far from the class center tend to have high leverage scores.

Additionally, we present a qualitative analysis of the images corresponding to the highest and lowest GLS values in Figure 3. We observe that high-leverage samples often exhibit atypical features or artifacts, such as unusual poses or backgrounds, which may contribute to their disproportionate influence on the model. Conversely, low-leverage samples tend to be more prototypical representations of their respective classes.

**Membership Inference Attack Evaluation** We evaluate the effectiveness of GLS in membership inference attacks by comparing the performance of LiRA for high-GLS samples and low-GLS samples. The results show that high-GLS samples are significantly more vulnerable to membership inference attacks compared to low-GLS samples.

Figure 4 illustrates the error trade-off curves for samples in the top and bottom 2% of GLS quantiles. The curves represent the mean of 50 LiRA attacks, with shaded areas indicating 95% percentile intervals. Notably, high-leverage samples demonstrate a marked shift toward the lower-left of the plot, indicating increased susceptibility to membership inference attacks.
<figure><img src="image_4.png" /><figcaption>Figure 3: Visualization of the 10 images with the Lowest and Highest Generalized Leverage Score.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Error trade-off curves ($\alpha$ = FPR, $\beta$ = FNR) for samples in the 2% highest and lowest GLS quantiles. Curves represent the mean of 50 LiRA attacks; shaded areas indicate 95% percentile intervals.</figcaption></figure>

# 7 Conclusion

We introduced the GLS, an influence-based, per-instance diagnostic for assessing privacy exposure without retraining models or explicitly simulating attacks. Although derived in the linear setting, GLS provides a computationally efficient tool for identifying data points with elevated risk of membership inference in deep learning models. Rather than offering formal privacy guarantees, it serves as an interpretable and scalable proxy for individual vulnerability that complements attack-based audits. While GLS does not perfectly align with shadow-model-based metrics, it consistently highlights samples that are more susceptible to privacy leakage, making it particularly suitable for large-scale auditing scenarios. A key limitation of our approach is that its theoretical justification is confined to linear models, with effectiveness in non-linear regimes supported empirically. Promising directions for future work include extending the framework to

white-box auditing settings, systematically evaluating its behavior across architectures and datasets, and studying its interaction with different training regimes, including differentially private optimization.

# References

* Martín Abadi, Andy Chu, Ian Goodfellow, H. Brendan McMahan, Ilya Mironov, Kunal Talwar, and Li Zhang. Deep Learning with Differential Privacy. In *Proceedings of the 2016 ACM SIGSAC Conference on Computer and Communications Security*, pages 308–318, October 2016. doi: 10.1145/2976749.2978318.

* David A Belsley, Edwin Kuh, and Roy E Welsch. *Regression Diagnostics: Identifying Influential Data and Sources of Collinearity*, chapter 2, pages 6–84. John Wiley & Sons, Ltd, 1980. doi: 10.1002/0471725153.ch2.

* Nicholas Carlini, Chang Liu, Úlfar Erlingsson, Jernej Kos, and Dawn Song. The secret sharer: Evaluating and testing unintended memorization in neural networks. In *28th USENIX security symposium (USENIX security 19)*, pages 267–284, 2019.

* Nicholas Carlini, Steve Chien, Milad Nasr, Shuang Song, Andreas Terzis, and Florian Tramer. Membership inference attacks from first principles. In *2022 IEEE symposium on security and privacy (SP)*, pages 1897–1914. IEEE, 2022.

* R. Dennis Cook. Detection of Influential Observation in Linear Regression. *Technometrics*, 19(1):15–18, 1977. doi: 10.1080/00401706.1977.10489493.

* Jinshuo Dong, Aaron Roth, and Weijie J Su. Gaussian differential privacy. *Journal of the Royal Statistical Society Series B: Statistical Methodology*, 84(1):3–37, 2022.
C. Dwork. *Differential privacy*, volume 2006. ICALP,
2006. Pages: 1-12.

Vitaly Feldman. Does learning require memorization?
a short tale about a long tail. In *Proceedings of the
52nd annual ACM SIGACT symposium on theory
of computing*, pages 954–959, 2020.

Vitaly Feldman and Chiyuan Zhang. What neural
networks memorize and why: Discovering the long
tail via influence estimation. *Advances in Neural In-
formation Processing Systems*, 33:2881–2891, 2020.

Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian
Sun. Deep residual learning for image recognition.
In *Proceedings of the IEEE conference on computer
vision and pattern recognition*, pages 770–778, 2016.

M.R. Hestenes and E. Stiefel. Methods of conjugate
gradients for solving linear systems. *Journal of
Research of the National Bureau of Standards*, 49
(6):409, December 1952. ISSN 0091-0635. doi: 10.
6028/jres.049.044.

Matthew Jagielski, Jonathan Ullman, and Alina
Oprea. Auditing differentially private machine learn-
ing: How private is private sgd? *Advances in Neural
Information Processing Systems*, 33:22205–22216,
2020.

Pang Wei Koh and Percy Liang. Understanding Black-
box Predictions via Influence Functions. In *Pro-
ceedings of the 34th International Conference on
Machine Learning*, pages 1885–1894. PMLR, July
2017.

Alex Krizhevsky and Geoffrey Hinton. Learning mul-
tiple layers of features from tiny images. Technical
Report 0, University of Toronto, Toronto, Ontario,
2009.

Mathias Lecuyer, Vaggelis Atlidakis, Roxana Geam-
basu, Daniel Hsu, and Suman Jana. Certified ro-
bustness to adversarial examples with differential
privacy. In *2019 IEEE symposium on security and
privacy (SP)*, pages 656–672. IEEE, 2019.

Milad Nasr, Shuang Songi, Abhradeep Thakurta, Nico-
las Papernot, and Nicholas Carlin. Adversary in-
stantiation: Lower bounds for differentially private
machine learning. In *2021 IEEE Symposium on
security and privacy (SP)*, pages 866–882. IEEE,
2021.

Jerzy Neyman and Egon Sharpe Pearson. Ix. on
the problem of the most efficient tests of statistical
hypotheses. *Philosophical Transactions of the Royal
Society of London, Series A: Containing Papers of a

*Mathematical or Physical Character*, 231(694-706):
289-337, 02 1933. doi: 10.1098/rsta.1933.0009.

Adam Paszke, Sam Gross, Francisco Massa, Adam
Lerer, James Bradbury, Gregory Chanan, Trevor
Killeen, Zeming Lin, Natalia Gimelshein, Luca
Antiga, et al. Pytorch: An imperative style, high-
performance deep learning library. *Advances in
neural information processing systems*, 32, 2019.

Barak A. Pearlmutter. Fast exact multiplication by
the hessian. *Neural Computation*, 6(1):147–160, 01
1994. doi: 10.1162/neco.1994.6.1.147.

Daryl Pregibon. Logistic regression diagnostics. *The Annals of Statistics*, 9(4):705–724, 1981.

Reza Shokri, Marco Stronati, Congzheng Song, and
Vitaly Shmatikov. Membership inference attacks
against machine learning models. In *2017 IEEE
symposium on security and privacy (SP)*, pages 3–
18. IEEE, 2017.

Laurens van der Maaten and Geoffrey Hinton. Vi-
sualizing Data using t-SNE. *Journal of Machine
Learning Research*, 9(86):2579–2605, 2008.

Samuel Yeom, Irene Giacomelli, Matt Fredrikson,
and Somesh Jha. Privacy risk in machine learn-
ing: Analyzing the connection to overfitting. In
*2018 IEEE 31st computer security foundations sym-
posium (CSF)*, pages 268–282. IEEE, 2018.

Sajjad Zarifzadeh, Philippe Liu, and Reza Shokri.
Low-Cost High-Power Membership Inference At-
tacks. In Ruslan Salakhutdinov, Zico Kolter,
Katherine Heller, Adrian Weller, Nuria Oliver,
Jonathan Scarlett, and Felix Berkenkamp, editors,
*Proceedings of the 41st International Conference on
Machine Learning*. PMLR, 21–27 Jul 2024.

Chiyuan Zhang, Samy Bengio, Moritz Hardt, Ben-
jamin Recht, and Oriol Vinyals. Understanding
deep learning requires rethinking generalization. In
*International Conference on Learning Representa-
tions*, 2017.
# A Leverage Scores and Membership Vulnerability: Proofs in the Linear Gaussian Case

In this appendix, we provide detailed proofs for the results presented in Section 3, which analyzes the relationship between leverage scores and membership inference vulnerability in the context of linear models.

## A.1 Setup and Notations

We recall here the setup for the linear regression model and the associated residual distributions under membership and non-membership hypotheses.

Consider a fixed design matrix $X \in \mathbb{R}^{n \times d}$, where each row $x_i^\top$ represents a data point. We examine the multivariate linear model:

$$Y = X\Theta^* + E,$$

where $Y \in \mathbb{R}^{n \times m}$ is the response matrix, $\Theta^* \in \mathbb{R}^{d \times m}$ is the true parameter matrix, and $E \in \mathbb{R}^{n \times m}$ represents noise. We assume the noise is centered, independent and identically distributed (i.i.d.) Gaussian, such that for any row $i$, $E[E_i.] = 0$ and $\text{Cov}(E_i.) = \sigma^2 I_m$.

The Ordinary Least Squares (OLS) estimator is $\hat{\Theta} = (X^\top X)^{-1} X^\top Y$. The fitted values are $\hat{Y} = HY$, where $H = X(X^\top X)^{-1}X^\top$ is the hat matrix. For a specific data point $i$, the residual vector $r_i \in \mathbb{R}^m$ is defined as:

$$r_i = y_i - \hat{\Theta}^\top x_i = \sum_{j=1}^{n} (\delta_{ij} - h_{ij}) y_j,$$

where $\delta_{ij}$ is the Kronecker delta. The diagonal elements of the hat matrix, $h_{ii} = x_i^\top (X^\top X)^{-1} x_i$, are the leverage scores. They satisfy $0 \le h_{ii} \le 1$ and $\sum_{i=1}^n h_{ii} = d$, serving as a measure of the geometric influence of $x_i$ on the model's predictions (Belsley et al., 1980).

We consider the task of Membership Inference Attack (MIA), framed as a hypothesis test on a residual $r_i$. We distinguish between two scenarios for a point $i$:

* $\mathcal{H}_1 = \mathcal{H}_{train}$ (Member): The point $i$ was included in the training set.

* $\mathcal{H}_0 = \mathcal{H}_{test}$ (Non-member): The point $i$ is a fresh test point with the same features $x_i$ but a fresh noise realization $\tilde{\varepsilon}_i$.

## A.2 Residual Distributions

**Proposition 3.1 (Residual Distributions).** Under the Gaussian noise assumption, the residuals for the $i$-th data point under the member and non-member hypotheses follow distinct multivariate normal distributions:

$$r_i | \mathcal{H}_{train} \sim \mathcal{N}(0, \sigma^2(1 - h_{ii})I_m), \quad (3)$$

$$r_i | \mathcal{H}_{test} \sim \mathcal{N}(0, \sigma^2(1 + h_{ii})I_m). \quad (4)$$

Consequently, the squared norms of the residuals follow scaled Chi-squared distributions:

$$\|r_i\|^2 | \mathcal{H}_{train} \sim \sigma^2(1 - h_{ii})\chi^2(m), \quad (5)$$

$$\|r_i\|^2 | \mathcal{H}_{test} \sim \sigma^2(1 + h_{ii})\chi^2(m). \quad (6)$$

*Proof.* Under the membership hypothesis $\mathcal{H}_1$, the residual for point $i$ is:

$$\begin{align*}
r_i &= y_i - \hat{\Theta}^\top x_i \\
&= (\Theta^{*\top} x_i + \varepsilon_i) - \sum_{j=1}^{n} h_{ij} (\Theta^{*\top} x_j + \varepsilon_j) \\
&= \varepsilon_i - \sum_{j=1}^{n} h_{ij} \varepsilon_j \\
&= (1 - h_{ii}) \varepsilon_i - \sum_{j \neq i} h_{ij} \varepsilon_j
\end{align*}$$
Since the noise terms $\epsilon_j$ are i.i.d. Gaussian with mean 0 and covariance $\sigma^2 I_m$, the residual $\mathbf{r}_i$ is also Gaussian with mean 0 and covariance:

$$
\begin{align*}
\mathrm{Cov}(\mathbf{r}_i) &= \mathbb{E}[\mathbf{r}_i \mathbf{r}_i^\top] \\
&= (1 - h_{ii})^2 \sigma^2 I_m + \sum_{j \neq i} h_{ij}^2 \sigma^2 I_m \\
&= \sigma^2 \left( (1 - h_{ii})^2 + \sum_{j \neq i} h_{ij}^2 \right) I_m \\
&= \sigma^2 (1 - h_{ii}) I_m
\end{align*}
$$

where we used the property of the hat matrix that $\sum_{j=1}^n h_{ij}^2 = h_{ii}$. Under the non-membership hypothesis $\mathcal{H}_0$, the residual for point $i$ is:

$$
\begin{align*}
\mathbf{r}_i &= \tilde{y}_i - \hat{\Theta}^\top \mathbf{x}_i \\
&= (\Theta^{*\top} \mathbf{x}_i + \tilde{\mathbf{\epsilon}}_i) - \sum_{j=1}^n h_{ij} (\Theta^{*\top} \mathbf{x}_j + \mathbf{\epsilon}_j) \\
&= \tilde{\mathbf{\epsilon}}_i - \sum_{j=1}^n h_{ij} \mathbf{\epsilon}_j \\
&= \tilde{\mathbf{\epsilon}}_i - h_{ii} \mathbf{\epsilon}_i - \sum_{j \neq i} h_{ij} \mathbf{\epsilon}_j
\end{align*}
$$

Similarly, the covariance of $\mathbf{r}_i$ under $\mathcal{H}_0$ is:

$$
\begin{align*}
\mathrm{Cov}(\mathbf{r}_i) &= \mathbb{E}[\mathbf{r}_i \mathbf{r}_i^\top] \\
&= \sigma^2 I_m + h_{ii}^2 \sigma^2 I_m + \sum_{j \neq i} h_{ij}^2 \sigma^2 I_m \\
&= \sigma^2 \left( 1 + h_{ii}^2 + \sum_{j \neq i} h_{ij}^2 \right) I_m \\
&= \sigma^2 (1 + h_{ii}) I_m
\end{align*}
$$

Thus, we have:

$$
\begin{align*}
\mathbf{r}_i | \mathcal{H}_{\text{train}} &\sim \mathcal{N}(0, \sigma^2(1 - h_{ii})I_m) \\
\mathbf{r}_i | \mathcal{H}_{\text{test}} &\sim \mathcal{N}(0, \sigma^2(1 + h_{ii})I_m)
\end{align*}
$$

□

### A.3 Optimal Detection via Likelihood Ratio

To distinguish between $\mathcal{H}_1$ and $\mathcal{H}_0$, the Neyman-Pearson lemma identifies the Likelihood Ratio (LR) as the most powerful test.

**Proposition 3.2 (Optimal MIA Test).** Let $S_i(\mathbf{r}_i)$ be the log-likelihood ratio statistic for the $i$-th sample. The most powerful test at level $\alpha$ rejects the non-member hypothesis $\mathcal{H}_{test}$ (predicts membership) if $S_i(\mathbf{r}_i) > \gamma_\alpha$, where the sufficient statistic is given by:

$$S_i(\mathbf{r}_i) = \frac{m}{2} \ln \left( \frac{1+h_{ii}}{1-h_{ii}} \right) - \frac{h_{ii}}{\sigma^2(1-h_{ii}^2)} \|\mathbf{r}_i\|^2. \quad (7)$$

*Proof.* The likelihood ratio for the residual $\mathbf{r}_i$ is given by:

$$\Lambda(\mathbf{r}_i) = \frac{p(\mathbf{r}_i | \mathcal{H}_{\text{train}})}{p(\mathbf{r}_i | \mathcal{H}_{\text{test}})}$$
where $p(\mathbf{r}_i|\mathcal{H}_{\text{train}})$ and $p(\mathbf{r}_i|\mathcal{H}_{\text{test}})$ are the probability density functions of the multivariate normal distributions under the two hypotheses. Using the results from Section A.2, we have:

$$ 
\begin{aligned}
p(\mathbf{r}_i|\mathcal{H}_{\text{train}}) &= \frac{1}{(2\pi)^{m/2}|\sigma^2(1-h_{ii})I_m|^{1/2}} \exp\left(-\frac{1}{2}\mathbf{r}_i^\top(\sigma^2(1-h_{ii})I_m)^{-1}\mathbf{r}_i\right) \\
&= \frac{1}{(2\pi)^{m/2}(\sigma^2(1-h_{ii}))^{m/2}} \exp\left(-\frac{1}{2\sigma^2(1-h_{ii})}\|\mathbf{r}_i\|^2\right)
\end{aligned}
$$

and

$$ 
\begin{aligned}
p(\mathbf{r}_i|\mathcal{H}_{\text{test}}) &= \frac{1}{(2\pi)^{m/2}|\sigma^2(1+h_{ii})I_m|^{1/2}} \exp\left(-\frac{1}{2}\mathbf{r}_i^\top(\sigma^2(1+h_{ii})I_m)^{-1}\mathbf{r}_i\right) \\
&= \frac{1}{(2\pi)^{m/2}(\sigma^2(1+h_{ii}))^{m/2}} \exp\left(-\frac{1}{2\sigma^2(1+h_{ii})}\|\mathbf{r}_i\|^2\right)
\end{aligned}
$$

Substituting these expressions into the likelihood ratio, we obtain:

$$ 
\begin{aligned}
\Lambda(\mathbf{r}_i) &= \frac{\frac{1}{(2\pi)^{m/2}(\sigma^2(1-h_{ii}))^{m/2}} \exp\left(-\frac{1}{2\sigma^2(1-h_{ii})}\|\mathbf{r}_i\|^2\right)}{\frac{1}{(2\pi)^{m/2}(\sigma^2(1+h_{ii}))^{m/2}} \exp\left(-\frac{1}{2\sigma^2(1+h_{ii})}\|\mathbf{r}_i\|^2\right)} \\
&= \left(\frac{1+h_{ii}}{1-h_{ii}}\right)^{m/2} \exp\left(\frac{1}{2\sigma^2}\left(\frac{1}{1+h_{ii}} - \frac{1}{1-h_{ii}}\right)\|\mathbf{r}_i\|^2\right) \\
&= \left(\frac{1+h_{ii}}{1-h_{ii}}\right)^{m/2} \exp\left(-\frac{h_{ii}}{\sigma^2(1-h_{ii}^2)}\|\mathbf{r}_i\|^2\right)
\end{aligned}
$$

Thus, the log-likelihood ratio is:

$$ \log \Lambda(\mathbf{r}_i) = \frac{m}{2} \log \left( \frac{1+h_{ii}}{1-h_{ii}} \right) - \frac{h_{ii}}{\sigma^2(1-h_{ii}^2)} \|\mathbf{r}_i\|^2 $$

<figure><img src="image_6.png" /></figure>

## A.4 Error Curve Derivation

We derive here the error trade-off curve for the optimal membership inference test based on the residual norm $\|\mathbf{r}_i\|^2$. We provide an illustration of these theoretical curves for different $h_{ii}$ values in Figure 5. Higher leverage scores correspond to curves closer to the origin, indicating increased vulnerability to membership inference attacks. The curves are non-symmetric due to the differing variances of the residual distributions under the two hypotheses, leading to easier detection of non-members than members.

**Proposition 3.3 (MIA Errors Curve).** For the point $i$, the trade-off between the false alarm rate $\alpha$ (Type-I error) and the missed detection rate $\beta$ (Type-II error) is described by the curve:

$$ \beta_i(\alpha_i) = 1 - F_m \left( \frac{1+h_{ii}}{1-h_{ii}} F_m^{-1}(\alpha_i) \right). \qquad (8) $$

where $F_m$ is the Cumulative Distribution Function of the $\chi^2(m)$ distribution.

*Proof.* The type I error (false positive rate, i.e., non-member misclassified as member) is

$$ \alpha(t) = \mathbb{P}_{\mathcal{H}_{\text{test}}} (\|\mathbf{r}_i\|^2 \le t) = F_{\chi_m^2} \left( \frac{t}{\sigma^2(1+h_{ii})} \right), \qquad (A.22) $$

where $F_{\chi_m^2}$ is the cumulative distribution function (CDF) of the $\chi^2$ distribution with $m$ degrees of freedom. The type II error (false negative rate, i.e., member misclassified as non-member) is

$$ \beta(t) = 1 - \mathbb{P}_{\mathcal{H}_{\text{train}}} (\|\mathbf{r}_i\|^2 \le t) = 1 - F_{\chi_m^2} \left( \frac{t}{\sigma^2(1-h_{ii})} \right). \qquad (A.23) $$
<figure><img src="image_7.png" /><figcaption>Figure 5: Theoretical error trade-off curves for different leverage scores $h_{ii}$ for $m = 1$.</figcaption></figure>

Eliminating $t$ by injecting Equation (A.22) into Equation (A.23), we obtain:

$$ \beta_i(\alpha_i) = 1 - F_m \left( \frac{1 + h_{ii}}{1 - h_{ii}} F_m^{-1}(\alpha) \right), $$

where $F_m$ is the CDF of the $\chi^2(m)$ distribution.

<figure><img src="image_8.png" /></figure>

# B Generalized Leverage Score: Proofs

In this section, we provide a formal derivation of the Generalized Leverage Score (GLS) for arbitrary differentiable models, following the notation and conventions of the main paper. We also derive the proofs for the second-order derivatives required in the closed-form expression of the GLS for two common loss functions: quadratic loss and cross entropy loss.

## B.1 Setup

We restate the problem setup for completeness. Consider a dataset $\mathcal{D} = \{(x_i, y_i)\}_{i=1}^n$ and a model $f_\theta: \mathcal{X} \to \mathcal{Y}$ parameterized by $\theta \in \Theta \subset \mathbb{R}^p$. The empirical risk is given by

$$ \mathcal{L}(\theta) = \frac{1}{n} \sum_{i=1}^{n} \ell(f_\theta(x_i), y_i), $$

where $\ell$ is a twice-differentiable loss function. The optimal parameters are $\hat{\theta} = \arg\min_{\theta \in \Theta} \mathcal{L}(\theta)$. We recall the definition of the Generalized Leverage Score.

**Definition 4.1 (Generalized Leverage Score).** The Generalized Leverage Score (GLS) of the $i$-th sample is defined as the infinitesimal sensitivity of the model's prediction at $x_i$ to its own observed label $y_i$:

$$ \text{GLS}_i = \frac{\partial f_{\hat{\theta}}(x_i)}{\partial y_i} \in \mathbb{R}^{m \times m}. \qquad (12) $$

## B.2 Derivation via Implicit Differentiation

We recall the closed-form expression for the GLS.

**Proposition 4.2 (Closed-form GLS).** Assuming the averaged loss function $\mathcal{L}$ is twice differentiable and its Hessian $H_{\hat{\theta}}$ is invertible, the GLS for a univariate output is given by:

$$ \text{GLS}_i = -\frac{1}{n} \| \nabla_{\theta} f_{\hat{\theta}}(x_i) \|_{H_{\hat{\theta}}^{-1}}^2 \left( \frac{\partial^2 \ell}{\partial y_i \partial f} \right), \qquad (13) $$
where $\|\mathbf{v}\|_A^2 = \mathbf{v}^\top A \mathbf{v}$ denotes the squared Matrix norm. For multivariate outputs, this generalizes to:

$$ \text{GLS}_i = -\frac{1}{n} \mathbf{J}_\theta f_{\hat{\theta}}(\mathbf{x}_i) \mathbf{H}_{\hat{\theta}}^{-1} \mathbf{J}_\theta f_{\hat{\theta}}(\mathbf{x}_i)^\top \left( \frac{\partial^2 \ell}{\partial \mathbf{y} \partial f} \right), \quad (14) $$

where $\mathbf{J}_\theta f_{\hat{\theta}}(\mathbf{x}_i) \in \mathbb{R}^{m \times p}$ is the Jacobian of the model outputs with respect to the parameters, $\mathbf{H}_{\hat{\theta}}$ is the Hessian of the loss function $\mathcal{L}$, i.e. $\mathbf{H}_{\hat{\theta}} = \nabla_{\hat{\theta}}^2 \mathcal{L}(\hat{\theta})$, and $\ell$ denotes the pointwise loss evaluated at $(f_{\hat{\theta}}(\mathbf{x}_i), \mathbf{y}_i)$.

*Proof.* Let $m \in \mathbb{N}^*$ be the output dimension of the model. We seek to compute the Jacobian $\text{GLS}_i = \frac{\partial f_{\hat{\theta}}(\mathbf{x}_i)}{\partial \mathbf{y}_i}$, accounting for the fact that $\hat{\theta}$ depends implicitly on $\mathbf{y}_i$ through the empirical risk minimization. Applying the chain rule, we have:

$$ \frac{\partial f_{\hat{\theta}}(\mathbf{x}_i)}{\partial \mathbf{y}_i} = \mathbf{J}_\theta f_{\hat{\theta}}(\mathbf{x}_i) \frac{\partial \hat{\theta}}{\partial \mathbf{y}_i}, \quad (B.24) $$

where $\mathbf{J}_\theta f_{\hat{\theta}}(\mathbf{x}_i) \in \mathbb{R}^{m \times p}$ is the Jacobian of $f_{\hat{\theta}}(\mathbf{x}_i)$ with respect to $\theta$, evaluated at $\hat{\theta}$.

We now need to compute $\frac{\partial \hat{\theta}}{\partial \mathbf{y}_i}$. We use implicit differentiation on the optimality condition for $\hat{\theta}$. This condition states that the gradient of the empirical risk vanishes at $\hat{\theta}$:

$$ \nabla_{\theta} \mathcal{L}(\hat{\theta}) = 0. $$

Since $\hat{\theta}$ is implicitly a function of $\mathbf{y}_i$, we treat the optimality condition as an equation involving both $\theta$ and $\mathbf{y}_i$. Specifically, we express the first-order condition as $\nabla_{\theta} \mathcal{L}(\hat{\theta}(\mathbf{y}_i), \mathbf{y}_i) = 0$, making explicit that $\hat{\theta}$ depends on $\mathbf{y}_i$ through the optimization, and that $\mathcal{L}$ also depends directly on $\mathbf{y}_i$ via the $i$-th term in the empirical risk.

$$ \begin{aligned} & \frac{d}{d\mathbf{y}_i} \nabla_{\theta} \mathcal{L}(\hat{\theta}(\mathbf{y}_i), \mathbf{y}_i) = 0 \\ & \nabla_{\theta}^2 \mathcal{L}(\hat{\theta}) \frac{\partial \hat{\theta}}{\partial \mathbf{y}_i} + \frac{\partial}{\partial \mathbf{y}_i} \nabla_{\theta} \mathcal{L}(\hat{\theta}, \mathbf{y}_i) = 0, \end{aligned} \quad (B.25) $$

where $\nabla_{\theta}^2 \mathcal{L}(\hat{\theta}) = \mathbf{H}_{\hat{\theta}}$ is the Hessian of the loss and where the partial derivative $\frac{\partial}{\partial \mathbf{y}_i} \nabla_{\theta} \mathcal{L}$ is taken with respect to $\mathbf{y}_i$ only via the $i$-th term in the empirical risk, holding $\hat{\theta}$ fixed. In the following we derive a closed-form for $\frac{\partial}{\partial \mathbf{y}_i} \nabla_{\theta} \mathcal{L}(\hat{\theta}, \mathbf{y}_i)$, and thus we drop the dependence of $\hat{\theta}$ on $\mathbf{y}_i$ in the notation for clarity.

Substituting the definition of the empirical risk:

$$ \begin{aligned} \frac{\partial}{\partial \mathbf{y}_i} \nabla_{\theta} \mathcal{L}(\hat{\theta}, \mathbf{y}_i) &= \frac{\partial}{\partial \mathbf{y}_i} \left( \frac{1}{n} \sum_{j=1}^{n} \nabla_{\theta} \ell(f_{\hat{\theta}}(\mathbf{x}_j), \mathbf{y}_j) \right) \\ &= \frac{1}{n} \sum_{j=1}^{n} \frac{\partial}{\partial \mathbf{y}_i} \nabla_{\theta} \ell(f_{\hat{\theta}}(\mathbf{x}_j), \mathbf{y}_j). \end{aligned} $$

Since the dataset samples are independent, the loss for sample $j$ (where $j \neq i$) does not depend explicitly on $\mathbf{y}_i$. Therefore, $\frac{\partial}{\partial \mathbf{y}_i} \nabla_{\theta} \ell_j = 0$ for all $j \neq i$. The summation collapses to the single $i$-th term:

$$ \begin{aligned} \frac{\partial}{\partial \mathbf{y}_i} \nabla_{\theta} \mathcal{L}(\hat{\theta}, \mathbf{y}_i) &= \frac{1}{n} \frac{\partial}{\partial \mathbf{y}_i} \nabla_{\theta} \ell(f_{\hat{\theta}}(\mathbf{x}_i), \mathbf{y}_i) \\ &= \frac{1}{n} \frac{\partial}{\partial \mathbf{y}_i} (\mathbf{J}_{\theta} f_{\hat{\theta}}(\mathbf{x}_i)^{\top} \nabla_{f} \ell(f_{\hat{\theta}}(\mathbf{x}_i), \mathbf{y}_i)) \\ &= \frac{1}{n} \mathbf{J}_{\theta} f_{\hat{\theta}}(\mathbf{x}_i)^{\top} \frac{\partial^2 \ell}{\partial \mathbf{y} \partial f}, \end{aligned} $$
as $J_{\theta}f_{\hat{\theta}}(\boldsymbol{x}_i)$ does not depend on $\boldsymbol{y}_i$. Substituting back into Equation (B.25):

$$ \frac{\partial \hat{\theta}}{\partial \boldsymbol{y}_i} = -\frac{1}{n} \boldsymbol{H}_{\hat{\theta}}^{-1} J_{\theta} f_{\hat{\theta}}(\boldsymbol{x}_i)^{\top} \frac{\partial^2 \ell}{\partial \boldsymbol{y} \partial f}. $$

and then substituting into the chain rule expression for GLS$_i$ (Equation (B.24)):

$$ \text{GLS}_i = -\frac{1}{n} J_{\theta} f_{\hat{\theta}}(\boldsymbol{x}_i) \boldsymbol{H}_{\hat{\theta}}^{-1} J_{\theta} f_{\hat{\theta}}(\boldsymbol{x}_i)^{\top} \left( \frac{\partial^2 \ell}{\partial f \partial \boldsymbol{y}} \right). $$

For $m=1$, this reduces to the scalar case with $J_{\theta}f_{\hat{\theta}}(\boldsymbol{x}_i) = \nabla_{\theta}f_{\hat{\theta}}(\boldsymbol{x}_i)^{\top}$ and $\frac{\partial^2 \ell}{\partial f \partial \boldsymbol{y}}$ being a scalar, recovering the following result:

$$ \begin{aligned} \text{GLS}_i &= -\frac{1}{n} \left( \frac{\partial^2 \ell}{\partial f \partial \boldsymbol{y}} \right) \nabla_{\theta} f_{\hat{\theta}}(\boldsymbol{x}_i)^{\top} \boldsymbol{H}_{\hat{\theta}}^{-1} \nabla_{\theta} f_{\hat{\theta}}(\boldsymbol{x}_i) \\ &= -\frac{1}{n} \left( \frac{\partial^2 \ell}{\partial f \partial \boldsymbol{y}} \right) \| \nabla_{\theta} f_{\hat{\theta}}(\boldsymbol{x}_i) \|_{\boldsymbol{H}_{\hat{\theta}}^{-1}}^2. \end{aligned} $$

□

## B.3 Loss Function Derivatives

We provide explicit derivations for the mixed second derivative $\frac{\partial^2 \ell}{\partial \boldsymbol{y} \partial f}$ for two common loss functions: quadratic loss and cross entropy loss with softmax outputs.

**Proposition 4.3 (Multivariate Regression).** For the quadratic loss with targets $\boldsymbol{y} \in \mathbb{R}^m$ and predictions $f \in \mathbb{R}^m$, the second derivative is:

$$ \frac{\partial^2 \ell}{\partial \boldsymbol{y} \partial f} = -2 \boldsymbol{I}_m. \qquad (15) $$

*Proof.* For the quadratic loss, we have:

$$ \ell(f(\boldsymbol{x}), \boldsymbol{y}) = \|f(\boldsymbol{x}) - \boldsymbol{y}\|_2^2. $$

The first derivative with respect to $f$ is:

$$ \frac{\partial \ell}{\partial f} = 2(f(\boldsymbol{x}) - \boldsymbol{y}). $$

The second derivative is:

$$ \begin{aligned} \frac{\partial^2 \ell}{\partial \boldsymbol{y} \partial f} &= \frac{\partial}{\partial \boldsymbol{y}}(2(f(\boldsymbol{x}) - \boldsymbol{y})) \\ &= -2\boldsymbol{I}_m, \end{aligned} $$

where $\boldsymbol{I}_m$ is the $m \times m$ identity matrix.

□

**Proposition 4.4 (Cross-Entropy Derivatives).** For the cross-entropy loss with one-hot labels $\boldsymbol{y}$ and logits $f$, the second derivative is:

$$ \frac{\partial^2 \ell}{\partial \boldsymbol{y} \partial f} = -\boldsymbol{I}_m. \qquad (17) $$

*Proof.* For the cross entropy loss with softmax outputs, we have:

$$ \begin{aligned} \ell(f(\boldsymbol{x}), \boldsymbol{y}) &= -\sum_{c=1}^{m} y_c \log \left( \frac{e^{f_c(\boldsymbol{x})}}{\sum_{j=1}^{m} e^{f_j(\boldsymbol{x})}} \right) \\ &= -\sum_{c=1}^{m} y_c f_c(\boldsymbol{x}) + \log \left( \sum_{j=1}^{m} e^{f_j(\boldsymbol{x})} \right). \end{aligned} $$
The first derivative with respect to f is:

$$ \frac{\partial \ell}{\partial f_c} = -y_c + \frac{e^{f_c(x)}}{\sum_{j=1}^{m} e^{f_j(x)}} = p_c - y_c, $$

where $p_c$ is the predicted probability for class c. The second derivative is:

$$ \frac{\partial^2 \ell}{\partial f_c \partial y_d} = \frac{\partial}{\partial y_d}(p_c - y_c) = -\delta_{cd}, $$

where $\delta_{cd}$ is the Kronecker delta, equal to 1 if c = d and 0 otherwise. Thus, the second derivative matrix is:

$$ \frac{\partial^2 \ell}{\partial y \partial f} = -I_m, $$

where $I_m$ is the m × m identity matrix.

<figure><img src="image_9.png" /></figure>

**Remark B.1.** For vector-valued outputs, GLS$_i$ is an m × m matrix. Scalar summaries such as the trace, Frobenius norm, or Spectral norm can be used to obtain a single leverage score per sample. As explained in the main part we decided to use the trace in our experiments as all these operations led to similar results.

## B.4 Proof of Proposition 4.5

**Proposition B.2 (Binary Logistic Leverage Score).** Under the cross entropy loss, for a binary logistic regression model, the GLS in the probability space reduces to:

$$ \frac{\partial \hat{p}_i}{\partial y_i} = \sigma'(\hat{y}_i) \times \text{GLS}_i = w_i \mathbf{x}_i^\top (\mathbf{X}^\top \mathbf{W} \mathbf{X})^{-1} \mathbf{x}_i, \quad (B.26) $$

with $w_i = \hat{p}_i(1 - \hat{p}_i) = \sigma'(\hat{y}_i)$ and $\mathbf{W} = \text{diag}(w_1, \dots, w_n)$.

*Proof.* We consider the binary logistic regression model where the prediction is $\hat{p}_i = \sigma(f_{\hat{\theta}}(\mathbf{x}_i))$ with $f_{\hat{\theta}}(\mathbf{x}_i) = \mathbf{x}_i^\top \hat{\theta}$. The individual cross entropy loss is given by $\ell(\hat{y}_i, y_i) = -y_i \log(\hat{p}_i) - (1 - y_i) \log(1 - \hat{p}_i)$.

### 1. Gradient of the model output

The gradient of the linear predictor $f_{\hat{\theta}}(\mathbf{x}_i)$ with respect to the parameters $\theta$ is:

$$ \nabla_{\theta} f_{\hat{\theta}}(\mathbf{x}_i) = \mathbf{x}_i $$

### 2. Second derivative of the loss

We first find the derivative of the loss with respect to the linear output $f$. Recall that $\frac{\partial \hat{p}_i}{\partial f} = \sigma'(f) = \hat{p}_i(1 - \hat{p}_i)$.

$$ \frac{\partial \ell}{\partial f} = \frac{\partial \ell}{\partial \hat{p}_i} \frac{\partial \hat{p}_i}{\partial f} = \left( -\frac{y_i}{\hat{p}_i} + \frac{1 - y_i}{1 - \hat{p}_i} \right) \hat{p}_i (1 - \hat{p}_i) = \hat{p}_i - y_i $$

Taking the mixed partial derivative with respect to the label $y_i$:

$$ \frac{\partial^2 \ell}{\partial y_i \partial f} = \frac{\partial}{\partial y_i} (\hat{p}_i - y_i) = -1 $$

### 3. Hessian of the total loss

The Hessian $\mathbf{H}_{\hat{\theta}}$ for the averaged cross entropy loss $\mathcal{L} = \frac{1}{n} \sum_j \ell_j$ is:
$$H_{\hat{\theta}} = \nabla_{\theta}^2 \mathcal{L} = \frac{1}{n} \sum_{j=1}^{n} \frac{\partial^2 \ell_j}{\partial f^2} x_j x_j^{\top}$$

Since $\frac{\partial^2 \ell_j}{\partial f^2} = \frac{\partial}{\partial f}(\hat{p}_j - y_j) = \hat{p}_j(1 - \hat{p}_j) = w_j$, we have:

$$H_{\hat{\theta}} = \frac{1}{n} X^{\top} W X$$

### 4. Computing the GLS$_i$ (Logit Space)

Substituting these into your general GLS formula:

$$\text{GLS}_i = -\frac{1}{n}(-1) \|x_i\|_{(\frac{1}{n}X^{\top}W X)^{-1}}^2 = x_i^{\top} (X^{\top}W X)^{-1} x_i$$

### 5. Mapping to the Probability Space

To find the sensitivity in the probability space $\frac{\partial \hat{p}_i}{\partial y_i}$, we apply the chain rule:

$$\frac{\partial \hat{p}_i}{\partial y_i} = \frac{\partial \hat{p}_i}{\partial f} \times \frac{\partial f}{\partial y_i} = \sigma'(\hat{y}_i) \times \text{GLS}_i$$

Substituting $w_i = \sigma'(\hat{y}_i) = \hat{p}_i(1 - \hat{p}_i)$, we obtain:

$$\frac{\partial \hat{p}_i}{\partial y_i} = w_i x_i^{\top} (X^{\top} W X)^{-1} x_i$$

This completes the proof.

<figure><img src="image_10.png" /></figure>

## C GLS for Last Linear Layer

In this section, we derive the exact closed-form solution for our score matrix when differentiation is restricted to the last linear layer. We denote the Generalized Leverage Score, computed with differentiation restricted to the last linear layer, as GLS$_i|_{\theta_{sub}}$. We also denote $d$ the dimension of the feature vector output by the penultimate layer, and we augment it with a bias term to form $\tilde{g}_i = [g(x_i)^{\top}, 1]^{\top} \in \mathbb{R}^{d+1}$.

We analyze two cases separately: the quadratic loss for regression and the cross entropy loss for classification. The separation is necessary because the quadratic loss treats each output dimension independently, whereas the cross entropy loss couples all classes through the softmax function, preventing decoupled computation. Remark C.3, a the end of the section, provide a complexity analysis for both cases.

### C.1 Quadratic Loss (Regression)

**Proposition C.1** (Last Linear Layer GLS for quadratic loss). Let $\tilde{g}_i \in \mathbb{R}^{d+1}$ be the feature vector (embeddings and bias) and $\tilde{G}$ the matrix stacking all vectors $\tilde{g}_i$. The last-layer GLS for sample $i$ under the quadratic loss is given by:

$$\text{GLS}_i|_{\theta_{sub}} = \tilde{g}_i^{\top} (\tilde{G}^{\top} \tilde{G})^{-1} \tilde{g}_i. \qquad (C.27)$$

*Proof.* We consider a neural network model where the last layer is a linear transformation.

#### 1. Notation and Jacobian Structure

Let the model output be $f(\boldsymbol{x}) = \boldsymbol{W}g(\boldsymbol{x}) + \boldsymbol{b} \in \mathbb{R}^m$. We denote the flattened parameter vector of this layer as $\theta_{sub} = \text{vec}([\boldsymbol{W}, \boldsymbol{b}]) \in \mathbb{R}^{m(d+1)}$. Since the model is linear with respect to $\theta_{sub}$, the Jacobian of the output $f(\boldsymbol{x}_i)$ with respect to the parameters is a tensor product:

$$J_i|_{\theta_{sub}} = \nabla_{\theta_{sub}} f(\boldsymbol{x}_i) = \boldsymbol{I}_m \otimes \tilde{g}_i^{\top} \in \mathbb{R}^{m \times m(d+1)} \qquad (C.28)$$

where $\boldsymbol{I}_m$ is the identity matrix of size $m$, and $\otimes$ denotes the Kronecker product.
## 2. The Hessian Matrix

The Hessian of the loss $\mathcal{L}$ with respect to $\theta_{sub}$ is defined as $H|_{\theta_{sub}} = \frac{1}{n} \sum_{j=1}^{n} J_{j}^{\top}|_{\theta_{sub}}(\nabla_{f}^{2} \ell_{j})J_{j}|_{\theta_{sub}}$. For the quadratic loss, we have $\nabla_{f}^{2} \ell = 2I_{m}$.

Thus, the Hessian becomes:

$$ 
\begin{aligned}
H|_{\theta_{sub}} &= \frac{2}{n} \sum_{j=1}^{n} (I_m \otimes \tilde{g}_j)(I_m \otimes \tilde{g}_j^{\top}) \\
&= I_m \otimes \left( \frac{2}{n} \sum_{j=1}^{n} \tilde{g}_j \tilde{g}_j^{\top} \right) \\
&= I_m \otimes \left( \frac{2}{n} \tilde{G}^{\top} \tilde{G} \right)
\end{aligned}
$$

The inverse is given by :

$$ H|_{\theta_{sub}}^{-1} = I_m \otimes \left( \frac{n}{2} (\tilde{G}^{\top} \tilde{G})^{-1} \right) \qquad (C.29) $$

## 3. Derivation of the GLS

The Generalized Leverage Score matrix is given by

$$ GLS_i = \frac{1}{n} \cdot J_i H^{-1} J_i^{\top} \cdot \left( -\frac{\partial^2 \ell}{\partial y \partial f} \right) $$

We have the equivalent formulation using the restricted Jacobian and Hessian:

$$ GLS_i|_{\theta_{sub}} = \frac{1}{n} \cdot J_i|_{\theta_{sub}} \cdot H|_{\theta_{sub}}^{-1} \cdot J_i^{\top}|_{\theta_{sub}} \cdot \left( -\frac{\partial^2 \ell}{\partial y \partial f} \right) \qquad (C.30) $$

The sensitivity term for quadratic loss is $-\frac{\partial^2 \ell}{\partial y \partial f} = -(2I_m) = 2I_m$.

Injecting the expressions derived in Equations (C.28) and (C.29) into Equation (C.30) we have :

$$ 
\begin{aligned}
GLS_i|_{\theta_{sub}} &= \frac{2}{n} (I_m \otimes \tilde{g}_i^{\top}) [I_m \otimes \frac{n}{2} (\tilde{G}^{\top} \tilde{G})^{-1}] (I_m \otimes \tilde{g}_i) \\
&= I_m \otimes (\tilde{g}_i^{\top} (\tilde{G}^{\top} \tilde{G})^{-1} \tilde{g}_i)
\end{aligned}
$$

This simplifies to the scalar leverage score multiplying the identity matrix.

<figure><img src="image_11.png" /></figure>

This formulation allows computing the GLS for all samples in $O(nd^2 + d^3)$ in time. The $O(nd^2)$ comes from computing $\tilde{G}^{\top} \tilde{G}$ and $O(d^3)$ from inverting the $d \times d$ matrix. This is more efficient than the solution proposed in Algorithm 2 for moderate $d$ and large $n$, as typically encountered in practice.

## C.2 Cross Entropy Loss (classification)

**Proposition C.2 (Last-Layer GLS for Cross Entropy).** Let $p_j = \text{softmax}(f(x_j)) \in \mathbb{R}^m$ be the probability vector for sample $j$, and let $S_j = \text{diag}(p_j) - p_j p_j^{\top}$ be the Jacobian of the loss w.r.t the logits. Let $\tilde{g}_i \in \mathbb{R}^{d+1}$ be the feature vector (embeddings and bias). The last-layer GLS for sample $i$ under the cross entropy loss is the $m \times m$ matrix where the entry $(u, v)$ is given by:

$$ (GLS_i|_{\theta_{sub}})_{uv} = \frac{1}{n} \tilde{g}_i^{\top} (H|_{\theta_{sub}})^{-1} [\tilde{g}_i]_{[u,v]} \qquad (C.31) $$

<footnote>This excludes the cost of the forward pass to obtain $\tilde{g}_i$ as these are typically available from the training or inference procedure. We refer only here to the computational overhead of the GLS computation itself.</footnote>
where $H|_{\theta_{sub}} \in \mathbb{R}^{m(d+1) \times m(d+1)}$ is the Hessian, structured as a block matrix where the $(k,l)$-th block (for $1 \le k, l \le m$) is given by:

$$H|_{\theta_{sub}[k,l]} = \frac{1}{n} \sum_{j=1}^{n} (S_j)_{kl} \tilde{g}_j \tilde{g}_j^{\top} \in \mathbb{R}^{(d+1) \times (d+1)}. \qquad (C.32)$$

*Proof.* We consider the same neural network model where the last layer is a linear transformation.

1. Notation and Jacobian Structure

As in the regression case, let $f(x) = Wg(x) + b \in \mathbb{R}^m$ and $\theta_{sub} = \text{vec}([W, b])$. The Jacobian of the output with respect to the parameters remains identical:

$$J_i|_{\theta_{sub}} = I_m \otimes \tilde{g}_i^{\top} \in \mathbb{R}^{m \times m(d+1)} \qquad (C.33)$$

2. The Hessian Matrix

The Hessian of the loss with respect to the parameters is defined as:

$$H|_{\theta_{sub}} = \frac{1}{n} \sum_{j=1}^{n} J_j^{\top}|_{\theta_{sub}} (\nabla_f^2 \ell_j) J_j|_{\theta_{sub}}$$

For the cross entropy loss, $\nabla_f^2 \ell_j = S_j$, which is the Hessian of the loss w.r.t the logits.

Substituting the Jacobian structure, we can explicitly compute the structure of $H|_{\theta_{sub}}$. The matrix consists of $m \times m$ blocks, where each block captures the interaction between the parameters of class $k$ and class $l$. The $(k, l)$-th block of size $(d+1) \times (d+1)$ is given by:

$$H|_{[k,l]} = \frac{1}{n} \sum_{j=1}^{n} \left( J_j^{\top}|_{\theta_{sub}} (S_j) J_j|_{\theta_{sub}} \right)_{[k,l]} \\ = \frac{1}{n} \sum_{j=1}^{n} (S_j)_{kl} \tilde{g}_j \tilde{g}_j^{\top}$$

Unlike the quadratic case, $S_j$ is dense, so $H|_{\theta_{sub}}$ is not block-diagonal (i.e., $H|_{[k,l]} \neq 0$ for $k \neq l$). Consequently, it must be inverted as a full matrix.

3. Derivation of the GLS

The Generalized Leverage Score matrix is given by:

$$\text{GLS}_i|_{\theta_{sub}} = \frac{1}{n} \cdot J_i|_{\theta_{sub}} \cdot H|_{\theta_{sub}}^{-1} \cdot J_i^{\top}|_{\theta_{sub}} \cdot \left( -\frac{\partial^2 \ell}{\partial y \partial f} \right) \qquad (C.34)$$

The sensitivity term for cross entropy is $-\frac{\partial^2 \ell}{\partial y \partial f} = I_m$.

Injecting Equation (C.33) into Equation (C.34), we have:

$$\text{GLS}_i|_{\theta_{sub}} = \frac{1}{n} (I_m \otimes \tilde{g}_i^{\top}) H|_{\theta_{sub}}^{-1} (I_m \otimes \tilde{g}_i)$$

To find the $(u, v)$-th scalar entry of this resulting $m \times m$ matrix, we observe that the Kronecker structure of the Jacobian effectively projects the $(u, v)$-th block of the inverse Hessian onto the feature space:

$$(\text{GLS}_i|_{\theta_{sub}})_{uv} = \frac{1}{n} \cdot \tilde{g}_i^{\top} (\mathbf{H}|_{\theta_{sub}}^{-1})_{[u,v]} \tilde{g}_i$$

□

*Remark C.3 (Computational Complexity Comparison).* The computational cost differs significantly between the two losses due to the structure of the Hessian.

For the quadratic loss, the Hessian consists of $m$ identical diagonal blocks, allowing us to compute and invert a single $(d+1) \times (d+1)$ matrix. This yields a total complexity of $OCHR^2 + d^3)$.

Conversely, for the cross entropy loss, the softmax function couples all $m$ classes, forcing the construction and inversion of the full $m(d+1) \times m(d+1)$ Hessian. This results in a significantly higher complexity of $OCHR^2 + m^3 d^3)$.