# Lazy Estimation of Variable Importance for Large Neural Networks

Yue Gao111Department of Statistics, University of Wisconsin Madison, Madison,
WI, USA. , Abby Stevens222Department of Statistics, University of Chicago,
Chicago, IL, USA. , Rebecca Willet††footnotemark:  333Department of Computer
Sciences, University of Chicago, Chicago, IL, USA. , Garvesh
Raskutti††footnotemark:

###### Abstract

As opaque predictive models increasingly impact many areas of modern life,
interest in quantifying the importance of a given input variable for making a
specific prediction has grown. Recently, there has been a proliferation of
model-agnostic methods to measure variable importance (VI) that analyze the
difference in predictive power between a full model trained on all variables
and a reduced model that excludes the variable(s) of interest. A bottleneck
common to these methods is the estimation of the reduced model for each
variable (or subset of variables), which is an expensive process that often
does not come with theoretical guarantees. In this work, we propose a fast and
flexible method for approximating the reduced model with important inferential
guarantees. We replace the need for fully retraining a wide neural network by
a linearization initialized at the full model parameters. By adding a ridge-
like penalty to make the problem convex, we prove that when the ridge penalty
parameter is sufficiently large, our method estimates the variable importance
measure with an error rate of $O(\frac{1}{\sqrt{n}})$ where $n$ is the number
of training samples. We also show that our estimator is asymptotically normal,
enabling us to provide confidence bounds for the VI estimates. We demonstrate
through simulations that our method is fast and accurate under several data-
generating regimes, and we demonstrate its real-world applicability on a
seasonal climate forecasting example.

##  1 Introduction

As predictive modeling becomes ubiquitous across a wide swath of application
areas, it is especially critical to understand which variables contribute most
to making a particular prediction. Black-box machine learning methods are
insufficient in the face of algorithmic decision-making about things like
sentencing, healthcare, and education, and working toward developing more
interpretable methods is becoming more and more relevant [RR19, GMR+18].

Traditional statistical tools based on parametric models (e.g. p-values,
ANOVA) for VI inference are dissatisfying due to restrictive assumptions often
violated in modern datasets. Non-parametric extensions thus have been explored
[DS95]. In recent decades, many VI methods designed for modern deep learning
models have been investigated; most of these methods are gradient-based and
depend on the structure and the weights of nodes in a given specific neural
network [SGK19, STY17, STK+17, BBM+15]. Few statistically rigorous properties
are provided for these methods, and the VI definition is always intimately
attached to the network itself, making it hard to interpret in a model-
agnostic setting.

In a model-agnostic setting, a natural definition of VI that is independent of
the estimation procedure is to measure the loss of predictive power when the
variables of interest are deleted. To estimate such model-agnostic VI,
_retraining_ is the most widely used type of method, which involves training
separate models on the reduced data with the variables of interest deleted and
assessing the predictive skill difference [WGSC21, LGR+18, SvdLP14].
Retraining often acts as the best benchmark to evaluate other VI estimation
methods [H+19] due to its accuracy, yet it is computationally infeasible in
high-dimensional settings. Other methods for VI estimation include knockoff
methods [BC18, CFJL17] and Floodgate [ZJ21], which require the co-variate
distribution to be known. An alternative approach is to use a _dropout_ -type
method [CRG17]. Dropout is best-suited for assessing how much a variable
affects a predictive model, as opposed to our goal of assessing how much a
variable affects the response. Despite the resulting issues with VI estimation
accuracy, it is still widely used in practice as a proxy for VI due to its
computational tractability.

In this work, we propose a computationally efficient variable importance
estimation procedure for model-agnostic and distribution-free settings with
theoretical guarantees that leverages a lazy retraining framework inspired by
[COB20]. The key idea is to train a new model on the transformed training
data, akin to retraining, but on a linearized version of the model centered
around model parameters learned from the original (unreduced) training data.
We perform ridge regression on this linearized model in the gradient feature
space, meaning that our lazy retraining procedure can be computed very
quickly. The resulting method, when applied to wide neural network models,
admits error bounds that show it is nearly as accurate as full retraining,
while computationally it is nearly as fast as dropout. Our theoretical bounds
are complemented by a collection of simulations that explore the limitations
of dropout and benefits of lazy retraining under a variety of conditions and
an application to understanding the importance of various climate indices in a
seasonal forecasting task.

In summary, the main contribution of this paper is a new, computationally
efficient VI estimation method with statistical performance guarantees in a
model-agnostic and distribution-free setting when using large neural networks.
Our theoretical analysis facilitates statistical inference, and we illustrate
our approach on both synthetic and real-world data to support the theoretical
claims and demonstrate the utility of our method. Other empirically-driven VI
estimation methods exhibit similarities to our approach; our theoretical
analyses may provide new insights into those methods as well as the one we
propose in this paper.

##  2 Notation and Preliminaries

Suppose we have samples $Z_{i}=({\mathbf{X}}_{i},Y_{i}),i=1,\dots,n$ for data
$Z=(X,Y)\sim P_{0}$, where ${\mathbf{X}}_{i}\in\mathbb{R}^{p}$ is the $i$-th
$p$-dimensional feature vector and $Y_{i}$ is the $i$-th observed response.
$X$ denotes the multi-variate random variable containing features, $Y$ denotes
the response random variable. Let $X_{-j}\in\mathbb{R}^{p-1}$ (resp.
${\mathbf{X}}_{i,-j}$) denote the features in $X$ (resp. ${\mathbf{X}}_{i}$)
with the $j$-th variable removed; on the other hand, if we replace the $j$-th
random variable in $X$ (resp. ${\mathbf{X}}_{i}$) by its marginal mean
$\mu_{j}=\mathbb{E}(X_{j})$, we denote it as $X^{(j)}$ (resp.
${\mathbf{X}}_{i}^{(j)}$), _i.e.,_
$X^{(j)}=(X_{1},\dots,X_{j-1},\mu_{j},X_{j+1},\dots,X_{p})$. Let
$P_{0},P_{0,-j}$ be the population distributions for $X$ and $X_{-j}$ and let
$P_{n},P_{n,-j}$ be the empirical distributions of $X$ and $X_{-j}$ for
$j\in[p]$. $\delta_{Z_{i}}$ denotes the point mass probability measure at the
$i$-th observation $Z_{i}$. We denote $\mathbb{E}_{0},\mathbb{E}_{0,-j}$ as
the expectations taken with respect to $P_{0}$ and $P_{0,-j}$.

Let $f_{0}$ denote the true function mapping $X$ to the expected value of $Y$
conditional on $X$, and let $f_{0,-j}$ denote the function mapping $X^{(j)}$
to the expected value of $Y$ conditional on $X^{(j)}$:

$\displaystyle f_{0}(X):=$ $\displaystyle\mathbb{E}_{0}[Y|X];$ (1)
$\displaystyle f_{0,-j}(X^{(j)}):=$
$\displaystyle\mathbb{E}_{0,-j}[Y|X_{-j}].$ (2)

Let $f_{n}$ be the empirical model trained using all $p$ variables in $X$
within a certain function class $\mathcal{F}$(we refer to this as the full
model):

$f_{n}\in\mathop{\rm
arg\min}_{f\in\mathcal{F}}\frac{1}{n}\sum_{i=1}^{n}[Y_{i}-f({\mathbf{X}}_{i})]^{2}.$
(3)

To measure the accuracy of an approximation $f_{n}(x)$ to its target function
$f_{0}$, we use the $L_{2}(\mu)$-norm

$\|f_{n}-f\|^{2}=\int|f_{n}(x)-f(x)|^{2}d\mu(x),$ (4)

where $\mu$ is the probability measure for $X$.

Further, we use $\epsilon$ and $\epsilon^{(j)}$ to denote the respective
remainder terms for any $j\in[p]$:

$\epsilon:=Y-\mathbb{E}_{0}[Y|X];\leavevmode\nobreak\ \leavevmode\nobreak\
\epsilon^{(j)}:=Y-\mathbb{E}_{0,-j}[Y|X_{-j}].$ (5)

We will define our measure of variable importance (VI) in terms of a
predictive skill measure $V(f,P)$ (the same measure in [WGSC21]). Larger
values of $V(f,P)$ should indicate better predictive performance. For
$Z=(X,Y)$, we denote $\dot{V}(f,P;\delta P)$ as the Gateaux derivative of
$V(f,P)$ at $P$ in the direction $\delta P$. Specifically, one of the
predictive skill measures we consider is the negative mean squared error
(MSE):

$V(f,P)=-\mathbb{E}_{(X,Y)\sim P}[Y-f(X)]^{2},$ (6)

and the corresponding $\dot{V}(f,P;\delta P)$ is $\dot{V}(f,P;\delta
P)=-\int_{Z=(X,Y)}(Y-f(X))^{2}d(\delta P)$. Hence, the Gateaux derivative of
the negative MSE is
$\dot{V}(f,P_{0};\delta_{Z_{i}}-P_{0})=-(Y_{i}-f(X_{i}))^{2}+\mathbb{E}[Y-f(X)]^{2}$
and $\mathbb{E}[\dot{V}(f_{0},P_{0};\delta_{Z_{i}}-P_{0})]=0$.

##  3 Estimating Variable Importance

The VI measure we consider, which makes no assumptions on the data generating
mechanism, is

${\text{\sc vi}_{j}}:=V(f_{0},P_{0})-V(f_{0,-j},P_{0,-j}).$ (7)

$\text{\sc vi}_{j}$ quantifies the difference in predictive skill between the
full model and the reduced model for any $j\in[p]$. Consider the following
simple linear model example, where we take the negative MSE as the
predictiveness measure.

###### Example 3.1

Suppose $Y=\beta_{1}X_{1}+\beta_{2}X_{2}+\epsilon$, where
$X_{i}\sim\mathcal{N}(0,\sigma^{2}),i=1,2,\leavevmode\nobreak\
\text{Cov}(X_{1},X_{2})=\rho$, and $\epsilon$ is a
$\mathcal{N}(0,\sigma_{\epsilon}^{2})$ noise that is independent of the
features. The variable importance of the first variable is

$\text{\sc
vi}_{1}=\beta_{1}^{2}\cdot\text{Var}(X_{1}|X_{2})=\beta_{1}^{2}(1-\rho^{2})\sigma^{2}$

due to the fact that $X_{1}|X_{2}\sim\mathcal{N}(\rho
X_{2},(1-\rho^{2})\sigma^{2})$ (see Section A.2.6)

In general, we see from this example that the variable importance measure is
determined not only by the relationship between $X_{j}$ and $Y$, but also the
covariance structure in the features.

Our goal is to estimate $\text{\sc vi}_{j}$ for any variable $X_{j}$ from data
$\\{({\mathbf{X}}_{i},Y_{i})\\}_{i=1}^{n}$ with no assumptions on the
relationship between $X$ and $Y$. For empirical estimators $f_{n}$ and
$f_{n,-j}$ of $f_{0}$ and $f_{0,-j}$, a plug-in estimator of our VI measure is

$\widehat{\text{\sc vi}}_{j}=V(f_{n},P_{n})-V(f_{n,-j},P_{n,-j}).$ (8)

The key problem we are concerned with in this paper is how to estimate
$f_{n,-j}$ in an accurate and computationally efficient way. Traditionally,
people use the following two types of methods to do the estimation: dropout
and retraining.

###  3.1 Dropout

The method we are calling dropout estimates $\mathbb{E}(Y|X_{-j})$ by plugging
the dropout features $X^{(j)}$ into the full model $f_{n}$. In this case, the
variable importance measure can be estimated by

$\widehat{\text{\sc vi}}^{(\text{\sc
dr})}_{j}=V(f_{n},P_{n})-V(f_{n},P_{n,-j}).$ (9)

For the negative MSE measure of predictive skill for instance, the dropout
estimate measures the difference between the squared error on the original
training set and the squared error on the training set after replacing feature
$j$ with its mean. Dropout is superior among all plug-in estimators in terms
of computational cost – we only need to train the model once to get $f_{n}$.
This is desirable, especially when the function class $\mathcal{F}$ is large
and complicated, such as with neural networks, and the computational cost for
training the model is high. Despite this benefit, dropout is unreliable in
many settings, as we will revisit in Section 3.3.

###  3.2 Retrain

An alternative to dropout is what we call retraining. Given a function class
$\mathcal{F}$, the retraining method estimates $\text{\sc vi}_{j}$ by training
separate models

${f}_{n,-j}\in\mathop{\rm
arg\min}_{f\in\mathcal{F}}[Y_{i}-f({\mathbf{X}}_{i}^{(j)})]^{2}$ (10)

for each variable $j\in[p]$ to estimate $f_{0,-j}$. Hence, VI under this
framework is estimated via

$\widehat{\text{\sc vi}}^{(\text{\sc
rt})}_{j}=V(f_{n},P_{n})-V(f_{n,-j},P_{n,-j}).$ (11)

When taking negative MSE as the predictive skill measure, the retraining
estimate in this case measures the difference between the squared error of a
model trained without feature $j$ and the squared error of a model trained
with feature $j$. Retraining is more accurate than dropout as long as the
function class $\mathcal{F}$ is large enough, but requires training $p+1$
models, which can be prohibitively computationally expensive in many settings.
In this paper, we are especially interested in the setting when the function
class is as large as a wide neural network.

###  3.3 Dropout vs. Retrain for Linear Models

The dropout method is widely used to estimate variable importance due to its
efficiency. However, in cases where variables in $X$ are highly correlated,
dropout behaves problematically. Below, we will illustrate and quantify the
difference of the variable importance estimation in the random design linear
model case, where we take the negative MSE as the $V(f,P)$ measure as in
Equation 6. For simplicity, we restrict the function space $\mathcal{F}$ to
the linear function space here.

Suppose $X\in\mathbb{R}^{p}\sim\mathcal{N}(0,\Sigma)$,
$\epsilon\sim\mathcal{N}(0,\sigma_{\epsilon}^{2})$. Assume $\Sigma$ is
positive definite. Let $\beta^{*}:=\mathop{\rm
arg\min}_{w\in\mathbb{R}^{p}}\mathbb{E}[Y-X^{\top}w]^{2}$, so
$\beta^{*}=\Sigma^{-1}\mathbb{E}(XY)$. In the population version, the dropout
method uses the predictor $X_{-j}^{\top}\beta^{*}_{-j}$ (where
$\beta^{*}_{-j}\in\mathbb{R}^{p-1}$ is $\beta^{*}$ with its $j$-th element
removed) to estimate $\mathbb{E}(Y|X_{-j})$, while the retraining method uses
the predictor $X_{-j}^{\top}{\beta}^{(j)}$, where
${\beta}^{(j)}\in\mathbb{R}^{p-1}$ is ${\beta}^{(j)}=\mathop{\rm
arg\min}_{w\in\mathbb{R}^{p-1}}\mathbb{E}[Y-X_{-j}^{\top}w]^{2}$. The
following proposition characterizes the difference between VI estimates
corresponding to the retraining and dropout methods.

###### Proposition 3.2

In the linear function space, the difference between the variable importance
estimates for variable $j$ from the population version of the dropout and
retraining methods is:

$\begin{split}&\widehat{\text{\sc vi}}^{(\text{\sc
dr})}_{j}-\widehat{\text{\sc vi}}^{(\text{\sc rt})}_{j}\\\
=&\frac{\vec{\bm{\gamma}_{j}}^{\top}\Sigma_{(j)}^{-1}\vec{\bm{\gamma}_{j}}}{(\Sigma_{jj}-\vec{\bm{\gamma}_{j}}^{\top}\Sigma_{(j)}^{-1}\vec{\bm{\gamma}_{j}})^{2}}\left[\mathbb{E}(X_{j}Y)-\vec{\bm{\gamma}_{j}}^{\top}\Sigma_{(j)}^{-1}\mathbb{E}(X_{-j}Y)\right]^{2},\end{split}$

where $\vec{\bm{\gamma}_{j}}=\mathbb{E}(X_{j}X_{-j})\in\mathbb{R}^{p-1}$.

If the true model between $Y$ and $X$ is linear, _i.e.,_
$Y=X^{\top}\beta^{*}+\epsilon$, and $X\perp\\!\\!\\!\perp\epsilon$, the
variable importance estimated by retraining linear regression is:

$\widehat{\text{\sc vi}}^{(\text{\sc
rt})}_{j}={\beta_{j}^{*}}^{2}(\Sigma_{jj}-\vec{\bm{\gamma}_{j}}^{\top}\Sigma_{(j)}^{-1}\vec{\bm{\gamma}_{j}});$
(12)

furthermore, in this setting $\widehat{\text{\sc vi}}^{(\text{\sc rt})}_{j}$
is exactly the true variable importance defined in (7). In contrast, the
dropout framework will give

$\widehat{\text{\sc vi}}^{(\text{\sc
dr})}_{j}={\beta_{j}^{*}}^{2}\cdot\Sigma_{jj}.$ (13)

If feature $j$ is important and highly correlated with feature $k$ (but
independent of all other features), then
$\vec{\bm{\gamma}_{j}}^{\top}\Sigma_{(j)}^{-1}\vec{\bm{\gamma}_{j}}$ may be
very large, making the difference between $\widehat{\text{\sc vi}}^{(\text{\sc
dr})}_{j}$ and $\widehat{\text{\sc vi}}^{(\text{\sc rt})}_{j}$ similarly
large. This example illustrates how dropout can significantly overestimate
variable importance, even in simple settings.

##  4 Lazy Training

Our central interest is in inferring VI using complex models that are time-
consuming to train, making the baseline retraining method described above
computationally infeasible. With this in mind, we turn our attention to neural
network (NN) models, a setting in which dropout is widely used.

Motivated by the need for faster and more accurate methods for estimating VI
with NN, we propose a computationally efficient VI estimate inspired by the
lazy training framework of [COB20] that estimates the difference between the
full model parameters and the model parameters when the $j$-th variable is
removed. Like dropout, our procedure only requires us to train the NN once on
the full data, and then we solve a linear system to update the full model
parameters for each variable $j\in[p]$.

Given the training data $\\{({\mathbf{X}}_{i}^{(j)},Y_{i})\\}$ sampled from
$(X,Y)\sim P_{0}$ for $i=1,\dots,n$ and the underlying function
$f_{0}(X)=\mathbb{E}_{P_{0}}[Y|X]$, there exists a a neural network function
class
$\\{h_{\theta}(x):\mathbb{R}^{p}\mapsto\mathbb{R}|\theta\in\mathbb{R}^{M}\\}$
that is parameterized by a vector $\theta$, such that when we train the model
parameters over this class by

$\theta_{f}=\mathop{\rm
arg\min}_{\theta\in\mathbb{R}^{M}}\frac{1}{n}\sum_{i=1}^{n}[Y_{i}-h_{\theta}({\mathbf{X}}_{i})]^{2},$
(14)

the estimation error can be bounded by
$\|h_{\theta_{f}}(x)-f_{0}(x)\|=O(n^{-1/2})$ up to some log terms [Bar94]. To
achieve this, the scale of the number of parameters $M$ depends on the
complexity of the target function.

For very complex functions, we can still achieve this accuracy with
$M=O(\sqrt{n})$.

In order to estimate $\text{\sc vi}_{j}$, we need an estimate of what we are
calling the reduced model $h_{\theta_{-j}}$, where

$\displaystyle\theta_{-j}=\mathop{\rm
arg\min}_{\theta\in\mathbb{R}^{M}}\frac{1}{n}\sum_{i=1}^{n}[Y_{i}-h_{\theta}({\mathbf{X}}_{i}^{(j)})]^{2}.$
(15)

Instead of retraining a NN to estimate $\theta_{-j}$, we can instead estimate
the difference between the full model parameters $\theta_{f}$ and
$\theta_{-j}$ using this linear approximation, and simply update the full
model parameters with this correction to estimate $h_{\theta_{-j}}$. We are
essentially regressing the error resulting from the dropout estimation against
the gradient to estimate this correction, and to do so we solve the following
convex problem based on the training data
$\\{({\mathbf{X}}_{i}^{(j)},Y_{i})\\}$ for $i=1,\dots,n$ and a $2$-norm
penalty on the parameters:

$\displaystyle\Delta\theta_{j}(\lambda,n)=\mathop{\rm
arg\min}_{\omega\in\mathbb{R}^{M}}$
$\displaystyle\Big{\\{}\frac{1}{n}\sum_{i=1}^{n}\big{[}Y_{i}-h_{\theta_{f}}({\mathbf{X}}_{i}^{(j)})$
(16)
$\displaystyle-\omega^{\top}\nabla_{\theta}h_{\theta}({\mathbf{X}}_{i}^{(j)})|_{\theta=\theta_{f}}\big{]}^{2}+\lambda\|\omega\|_{2}^{2}\Big{\\}},$

where $\lambda>0$ is the penalty parameter.

Accordingly, the reduced neural network parameters are
$\Delta\theta_{j}(\lambda,n)+\theta_{f}$. For the simplicity of notation, we
write $\Delta\theta_{j}(\lambda,n)$ as $\Delta\theta_{j}$ for short. Then the
reduced model approximation without the $j$-th feature is
$\mathbb{R}^{p}\mapsto\mathbb{R}:x\mapsto h_{\theta_{f}+\Delta\theta_{j}}(x)$.

Hence, the variable importance measure under lazy training is

$\widehat{\text{\sc vi}}^{(\text{\sc
lazy})}_{j}=V(h_{\theta_{f}},P_{n})-V(h_{\theta_{f}+\Delta\theta_{j}},P_{n,-j}).$
(17)

Under the negative MSE measure $V(f,P)$, we have

$\widehat{\text{\sc vi}}^{(\text{\sc
lazy})}_{j}=\frac{1}{n}\sum_{i=1}^{n}\\{[Y_{i}-h_{\theta_{f}+\Delta\theta_{j}}({\mathbf{X}}_{i}^{(j)})]^{2}-[Y_{i}-h_{\theta_{f}}({\mathbf{X}}_{i})]^{2}\\}.$

(More precisely, we use data splitting for training and estimating VI as
detailed in Algorithm 1.) Essentially, the linearized approximation of the NN
is linear in the gradient feature map
$x\mapsto\nabla_{\theta}h_{\theta}(x)|_{\theta_{f}}$. In fact, this gradient
feature map induces the Neural Tangent Kernel (NTK, [JGH20]): for any
$x,x^{\prime}\in\mathbb{R}^{p}$,

$\text{ker}_{\theta_{f}}(x,x^{\prime}):=\langle\nabla_{\theta}h_{\theta}(x)|_{\theta_{f}},\nabla_{\theta}h_{\theta}(x^{\prime})|_{\theta_{f}}\rangle.$
(18)

Thus $\Delta\theta_{j}$ can be viewed as the solution for a kernel ridge
regression problem with kernel $\text{ker}_{\theta_{f}}$.

###  4.1 Theoretical Guarantee

By [WGSC21], when the empirical estimates for $\mathbb{E}(Y|{\mathbf{X}})$ and
$\mathbb{E}(Y|{\mathbf{X}}_{-j})$ converge to the target functions $f_{0}$ and
$f_{0,-j}$ at the rate of $O_{p}(n^{-1/4})$ in function norm, we achieve an
asymptotically normal and efficient estimator for the VI measure. In this
section, we give a theoretical guarantee to show that the lazy prediction
$h_{\theta_{f}+\Delta\theta_{j}}(X^{(j)})$ for the reduced model achieves such
convergence rate, so that the lazy training procedure gives an accurate
estimate of VI with an error in the order of $O(\frac{1}{\sqrt{n}})$ and we
can make inference accordingly.

Let ${\mathbf{e}}^{(j)}$ denote the difference between the true reduced
function $f_{0,-j}(X^{(j)})$ and the corresponding dropout estimation:

${\mathbf{e}}^{(j)}:=f_{0,-j}({\mathbf{X}}^{(j)})-h_{\theta_{f}}({\mathbf{X}}^{(j)})\in\mathbb{R}^{n}.$
(19)

Further, we denote the kernel matrix on $X^{(j)}$ induced by the gradient
feature map as ${\mathbb{K}}^{(j)}\in\mathbb{R}^{n\times n}$, whose elements
are defined as:

${\mathbb{K}}^{(j)}_{ik}:=\text{ker}_{\theta_{f}}({\mathbf{X}}_{i}^{(j)},{\mathbf{X}}_{k}^{(j)}),\leavevmode\nobreak\
i,k\in[n].$ (20)

Before diving into the main results, we first clarify two types of notation
for order of approximation:

  * •

$f(n)=O(g(n))$ if there exists $M>0$ and $N>0$, such that $|f(n)|\leq Mg(n)$
for all $n>N$.

  * •

$X_{n}=O_{p}(a_{n})$ as $n\rightarrow\infty$ if for any $\epsilon>0$, there
exists $M>0$ and $N>0$, such that
$\mathbb{P}(|\frac{X_{n}}{a_{n}}|>M)<\epsilon$ for any $n>N$.

###### Assumption 4.1

For any $j\in[p]$ and the regularization parameter $\lambda=O(\sqrt{n})$, we
assume:

  * (a)

$\|[{\mathbb{K}}^{(j)}+\lambda
I_{n}]^{-1}{\mathbf{e}}^{(j)}\|^{2}=O_{p}(1/\sqrt{n})$;

  * (b)

$\text{tr}({\mathbb{K}}^{(j)})=O_{p}(n)$.

The above assumption (b) is commonly used in NTK literature (see e.g.
[HLY19]). For a two-layer neural network, we can verify this numerically (see
Section B.1). For the assumption (a), by the fact that ${\mathbb{K}}^{(j)}$ is
positive semi-definite, this assumption can be satisfied when we have a large
regularization $\lambda=O(\sqrt{n})$.

###### Assumption 4.2

For the noise term $\epsilon^{(j)}$, we have the following assumption on its
conditional tail probability: there exists $\sigma$ such that for any
$j\in[p]$,

$\mathbb{E}\left[e^{\lambda\epsilon^{(j)}}|X^{(j)}\right]\leq
e^{\sigma^{2}\lambda^{2}/2},\leavevmode\nobreak\ \text{ for all
}\lambda\in\mathbb{R}.$ (21)

###### Assumption 4.3

Denote the gradient feature matrix as $\Phi\in\mathbb{R}^{n\times
M}=(\nabla_{\theta}h_{\theta}({\mathbf{X}}_{1})|_{\theta=\theta_{f}},\dots,\nabla_{\theta}h_{\theta}({\mathbf{X}}_{n})|_{\theta=\theta_{f}})^{\top}$.
We assume $\|\Phi^{\top}{\mathbf{e}}^{(j)}\|_{2}\leq O_{p}(1)$.

This assumption essentially requires that the linear space of neural tangent
kernels can well represent ${\mathbf{e}}^{(j)}$. We know that
${\mathbf{e}}_{i}^{(j)}$ is a function of ${\mathbf{X}}_{i}^{(i)}$, thus as
long as the neural network function class is large enough, this can be
satisfied with respect to the sample size $n$.

###### Theorem 4.4

Suppose 4.1, 4.2 and 4.3 hold, then for a neural network structure
$h_{\theta}(\cdot)$ which is $L$-smooth with respect to its parameters
$\theta$, as long as we take the ridge penalty parameter in the order
$\lambda=O(n^{1/2})$, then the lazy training method can accurately predict the
reduced model without the $j$-th covariate, _i.e.,_

$\|h_{\theta_{f}+\Delta\theta_{j}}(x)-\mathbb{E}(Y|X^{(j)})\|_{2}=O_{p}(n^{-1/4}).$
(22)

Therefore our variable importance estimator $\widehat{\text{\sc
vi}}^{(\text{\sc lazy})}_{j}$ is asymptotically normal and has an error rate
$O_{p}(n^{-1/2})$:

$\widehat{\text{\sc vi}}^{(\text{\sc lazy})}_{j}-\text{\sc
vi}_{j}=\Delta_{n,j}+O_{p}(n^{-1/2});$ (23)

where

$\displaystyle\Delta_{n,j}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\big{[}\dot{V}(f_{0},P_{0};\delta_{Z_{i}}-P_{0})$
(24)
$\displaystyle-\dot{V}(f_{0,-j},P_{0,-j};\delta_{Z_{i}}-P_{0,-j})\big{]}\rightarrow_{d}\mathcal{N}(0,\tau_{n,j}^{2});$

here the variance is
$\tau_{n,j}^{2}=\text{Var}({\epsilon^{(j)}}^{2}-\epsilon^{2})/n$, where
$\epsilon$ and $\epsilon^{(j)}$ is defined in Equation 5.

This result enables us to construct Wald-type confidence intervals around our
LazyVI estimates. In particular, the $\alpha-$level confidence intervals are
given by

$\widehat{\text{\sc vi}}^{(\text{\sc lazy})}_{j}\pm
z_{\frac{\alpha}{2}}\hat{\tau}_{n,j}$ (25)

where $\hat{\tau}_{n,j}$ is the plug-in estimate of $\tau_{n,j}$ in (24) and
$z_{\frac{\alpha}{2}}$ is the $\alpha/2$ quantile of the standard normal
distribution.

###  4.2 Proof Overview

The challenge of proving Theorem 4.4 is to bound the error of the lazy neural
network trained using data without a certain variable – note that we are
bounding the estimation error ($\|h_{\theta_{f}+\Delta\theta_{j}}-f_{0,-j}\|$)
instead of the prediction error ($\|h_{\theta_{f}+\Delta\theta_{j}}-f_{0}\|$)
that is the focus of much of the deep learning community, since the predictive
skill of the reduced model is expected to decrease when an important variable
is removed. At a high level, our proof reduces the estimation error of the
neural network from lazy training to the error between the NTK estimation and
the target function, where we use techniques from kernel ridge regression. The
difference here is that most NTK papers (see e.g. [JGH20]) use random
initialization for the parameters and optimization without penalty, while our
method starts from a specific initialization (the full model), and requires
the penalty parameter $\lambda$ to be large ($\lambda=O(n^{1/2})$) to ensure
convergence.

The following two lemmas give some intuition on how the neural network trained
by the lazy procedure can accurately estimate the reduced model. Basically,
the bound for the error consists of two parts: the error from the kernel ridge
regression (discussed in Lemma 4.5), and the error from the linear
approximation of the neural network (in Lemma 4.6). More proof details are
deferred to the Appendix.

Denote the linear approximation of the network as

$\tilde{h}_{\theta_{f}+\Delta\theta_{j}}(x):=h_{\theta_{f}}(x)+\langle\nabla_{\theta}h_{\theta}(x)|_{\theta=\theta_{f}},\Delta\theta_{j}\rangle.$
(26)

###### Lemma 4.5

Let $\lambda$ be penalty parameter in Equation 16, we have with probability at
least $1-\delta$,

$\displaystyle\|\tilde{h}_{\theta_{f}+\Delta\theta_{j}}(X^{(j)})-f_{0,-j}(X^{(j)})\|_{n}$
(27) $\displaystyle\leq\frac{\lambda\|[{\mathbb{K}}^{(j)}+\lambda
I_{n}]^{-1}{\mathbf{e}}^{(j)}\|}{\sqrt{n}}+\sigma\sqrt{\frac{\text{tr}[{\mathbb{K}}^{(j)}]}{4\lambda
n}}+\sigma\sqrt{\frac{2\log(1/\delta)}{n}}.$

Lemma 4.5 combined with 4.1 when the penalty parameter is
$\lambda=O(n^{1/2})$, yields a bound on the empirical error of the kernel
ridge regression component of $O_{p}(n^{-1/4})$. Based on this empirical
bound, we could then further bound the generalization error of the estimated
function using function complexity (See Section A.2.3).

###### Lemma 4.6

For a large neural network with width in the order $O(\sqrt{n})$, with high
probability we have for all $j\in[p]$,

$\|\tilde{h}_{\theta_{f}+\Delta\theta_{j}}(x)-h_{\theta_{f}+\Delta\theta_{j}}(x)\|\leq
O(n^{-1/4}).$ (28)

Lemma 4.6 shows that as long as the neural network is sufficiently large, the
neural network with updated parameters $\theta_{f}+\Delta_{j}$ is close to its
linear approximation.

###  4.3 Implementation

We estimate $h_{\theta_{f}}$ and $h_{\theta_{-j}}$ using $n_{1}<n$ samples as
training data, and use the remaining $n_{2}=n-n_{1}$ samples to estimate VI.
For the dropout method, VI is estimated simply by plugging the modified
testing data $\\{{\mathbf{X}}_{i}^{(j)}\\}_{i=n_{1}+1}^{n}$ into
$h_{\theta_{f}}$. For the retraining method, first $h_{\theta_{-j}}$ is
estimated by retraining the NN $h$ with
$\\{{\mathbf{X}}_{i}^{(j)}\\}_{i=1}^{n_{1}}$, and then VI is estimated by
plugging the modified testing data into this retrained estimate.

For the lazy training method, which we call LazyVI, we use the training data
to estimate the full model parameters, compute the gradient of the network
with respect to each model parameter for each training sample, and then
regress these gradients against the difference between $Y$ the dropout
estimates from the training data to estimate the parameter correction
$\Delta\theta_{j}$ for variable $j$. We then update the full model parameters
using this learned correction to compute the VI estimate and its associated
standard errors. See Algorithm 1 for full details.

Theorem 4.4 makes the assumption that the ridge parameter $\lambda$ from
Equation (16) is large. Since we are ultimately interested in estimating
$h_{\theta_{-j}}$ and not $\Delta\theta_{j}$, we evaluate
$h_{\theta_{f}+\Delta\theta_{j}}(\cdot)$ through K-fold CV to choose
$\hat{\lambda}_{j}$ for each variable (Algorithm 2 in Section B.2). Our
implementation is available at https://github.com/Willett-Group/lazyvi.

Algorithm 1 Lazy training for VI

0: Data: $\\{{\mathbf{X}}_{i},Y_{i}\\}_{i=1}^{n}$; $\lambda>0$; training size:
$0<n_{1}<n$; $n_{2}\leftarrow n-n_{1}$; NN structure:
$\theta\in\mathbb{R}^{M}\mapsto h_{\theta}(\cdot)$
LazyVI$\\{{\mathbf{X}}_{i},Y_{i}\\}_{i=1}^{n}$; $\lambda$, $n_{1}$

$\theta_{f}\leftarrow\mathop{\rm
arg\min}_{\theta\in\mathbb{R}^{M}}\frac{1}{n_{1}}\sum_{i=1}^{n_{1}}[Y_{i}-h_{\theta}({\mathbf{X}}_{i})]^{2}$

$v_{n}\leftarrow-\frac{1}{n_{2}}\sum_{i=n_{1}+1}^{n}[Y_{i}-h_{\theta_{f}}({\mathbf{X}}_{i})]^{2}$

for $j\in[p]$ do

${\mathbf{X}}_{i}^{(j)}\leftarrow{\mathbf{X}}_{i}$;
${\mathbf{X}}_{ij}^{(j)}\leftarrow\frac{1}{n_{1}}\sum_{i=1}^{n_{1}}{\mathbf{X}}_{ij}$

${\mathbf{e}}_{i}^{(j)}\leftarrow
Y_{i}-h_{\theta_{f}}({\mathbf{X}}_{i}^{(j)}),\leavevmode\nobreak\
i=1,\dots,n_{1}$

$\Phi_{i}^{(j)}\leftarrow\nabla_{\theta}h_{\theta}({\mathbf{X}}_{i}^{(j)})|_{\theta=\theta_{f}},\leavevmode\nobreak\
i=1,\dots,n_{1}$

$\Delta\theta_{j}\leftarrow\displaystyle{\mathop{\rm
arg\min}_{\omega\in\mathbb{R}^{M}}}\frac{1}{n_{1}}\sum_{i=1}^{n_{1}}[{\mathbf{e}}_{i}^{(j)}-\omega^{\top}\Phi_{i}^{(j)}]^{2}+\lambda\|\omega\|_{2}^{2}$

$v_{n,-j}\leftarrow\displaystyle-\frac{1}{n_{2}}\sum_{i=n_{1}+1}^{n}[Y_{i}-h_{\theta_{f}+\Delta\theta_{j}}({\mathbf{X}}_{i}^{(j)})]^{2}$

$\widehat{\text{\sc vi}}_{j}\leftarrow v_{n}-v_{n,-j}$

$t_{i,j}\leftarrow(Y_{i}-h_{\theta_{f}+\Delta\theta_{j}}({\mathbf{X}}_{i}^{(j)}))^{2}-(Y_{i}-h_{\theta_{f}}({\mathbf{X}}_{i}))^{2}$

$\hat{\tau}_{j}\leftarrow\frac{1}{n_{2}}\sum_{i=1}^{n_{2}}(t_{i,j}-\bar{t}_{j})^{2}/{n_{2}}$

end for

$\widehat{\text{\sc vi}}_{j},\leavevmode\nobreak\ j=1,\dots,p$.

##  5 Simulations

We first assess the performance of LazyVI on simulated data to highlight key
theoretical claims and assumptions and show that our method is empirically
practical. For these experiments, we train a wide, fully connected two-layer
neural network with ReLU activation for all simulations. Unless otherwise
specified, the width of the hidden layer in the training network is $m=50$.

###  5.1 Impact of Correlation in Linear Systems

Our first set of simulations serve to support key details of our theoretical
analysis. We consider data generated from the linear model
$f(X)=1.5X_{1}+1.2X_{2}+X_{3}+\epsilon$, where
$\epsilon\sim\mathcal{N}(0,0.1)$ and $X\sim\mathcal{N}(0,\Sigma_{6\times 6})$,
so the response only depends on the first three of the six variables. All
variables are independent except for $X_{1}$ and $X_{2}$, whose correlation is
$\rho$. As discussed in Example 3.1, the true VI of $X_{1}$, $X_{2}$, and
$X_{3}$ are given by $(1.5)^{2}(1-\rho^{2})$, $(1.2)^{2}(1-\rho^{2})$, and
$1$, respectively, and the VI of the remaining 3 variables is zero. In this
simple setting, we find that LazyVI approximates the true vi well with
desirable coverage and a considerable speed-up relative to retraining (Section
B.3).

We show in Prop. 3.2 that, when data are generated from a linear model, the
difference between the dropout and retraining variable importance estimates is
a function of the covariance of $X$. After training the full model, we use
both the dropout and our lazy procedure to estimate VI for increasing values
of $\rho$. In Figure 1, we show the difference between the dropout and LazyVI
estimates for variables $X_{1}$ and $X_{2}$ alongside the analytic difference
between $\widehat{\text{\sc vi}}^{(\text{\sc dr})}$ and vi (dotted line). We
see that the gap between LazyVI and dropout evolves with $\rho$ according to
the theoretical analysis, providing evidence that LazyVI behaves as expected.

Figure 1:  Difference between the dropout and LazyVI estimates for $X_{1}$ and
$X_{2}$. Dotted line is theoretical gap and shading shows std. across 10
repetitions. We see that LazyVI closely approximates the theoretical variable
importance measures under increasing correlations.

We use this simple linear setting to explore two additional assumptions from
our theoretical results. First, the linearization in (26) is a first order
Taylor approximation and assumes the full model parameters are close to the
reduced model parameters. If we try to linearize a neural network around a
random initialization, our LazyVI estimates are much less accurate and more
highly variable (Section B.4). Next, our theory assumes that our training
network is over-parameterized and sufficiently wide. We compute empirical
confidence intervals for LazyVI for increasing network widths and find that
coverage increases as the width increases, but at a computational cost
(Section B.5).

###  5.2 Binary Classification

Because we borrow much of our theoretical framework from [WGSC21], we also
leverage their simulation framework as a useful point of comparison. We draw
independent samples $X\sim\mathcal{N}(0,I_{4\times 4})$ and generate a binary
outcome $Y\sim\text{Bernoulli}(\Phi(X\beta))$ where $\beta=(2.5,3.5,0,0)$.
Because the outcome is binary, we use accuracy as our predictive skill
measure, and the true VI values are given by (0.136, 0.236, 0, 0),
respectively. We first directly compare the LazyVI and retrain estimators by
estimating vi across 100 simulated datasets of sample size $n=1000$ and
computing the empirical 95% confidence intervals. In Figure 2, we see that the
LazyVI and retrain estimates both achieve the desired level of coverage with
low bias. In this simulation, LazyVI took on average 0.6 seconds (including
cross-validating to find the optimal ridge parameter), while retraining took
$7.5$ seconds. In this setting, LazyVI is just as accurate as retraining with
a more than 10x speed-up.

Figure 2:  Left: Average coverage of empirical 95% confidence intervals from
the LazyVI and retrain estimates across 100 simulations. Right: Average
empirical bias ($\text{\sc vi}-\hat{\text{\sc vi}}$) of LazyVI and retrain
estimates. We see that LazyVI achieves low bias and high coverage in this
simulation setting, similar to the Retrain estimates.

###  5.3 Nonlinear, High-dimensional Regression

The computational burden of retraining is most pronounced in high-dimensional
settings, since estimating $\widehat{\text{\sc vi}}^{(\text{\sc rt})}$ for all
variables requires refitting at least $p$ models. For this simulation, we have
data $X\sim N(0,\Sigma_{100\times 100})$, where variables are independent
except $\text{Corr}(X_{1},X_{2})=0.5$. Letting
$\beta=(5,4,3,2,1,0,\dots,0)^{\top}\in\mathbb{R}^{100}$, we construct a weight
matrix $W\in\mathbb{R}^{m\times p}$ such that the
$W_{:,j}\sim\mathcal{N}(\beta_{j},\sigma^{2})$ (i.e. the weights associated
with variable $j$ are centered at $\beta_{j}$). Letting
$V\sim\mathcal{N}(0,1)$, we generate the response
$Y_{i}=V\sigma(W{\mathbf{X}}_{i})+\epsilon_{i}$ where $\sigma$ is the ReLU
function. Because the “true” VI values are unknown and difficult to estimate,
we present the accuracy of different estimation methods relative to the
retraining estimates, which we take as ground truth. We estimate vi for
$X_{1}$ across 10 simulated datasets ($n=1000$) and benchmark against
retraining using both a linear regression (OLS) and random forest (RF).In
Figure 3, we show the spread of both the computation time and normalized error
(relative to retrain) for all methods. We see that LazyVI is the most accurate
method and is substantially faster than retraining, which is especially
beneficial in this high-dimensional setting.

Figure 3: Distribution of computation time vs. normalized estimation error
relative to retrain for the VI of $X_{1}$ (($\hat{\text{\sc
vi}}-\widehat{\text{\sc vi}}^{(\text{\sc rt})})/\widehat{\text{\sc
vi}}^{(\text{\sc rt})})$ across 10 repetitions. Among the non-Retrain methods,
LazyVI is both the fastest and most accurate.

##  6 Predicting Seasonal Precipitation

Extreme precipitation events have become more and more common in recent years,
and are expected to intensify with climate change [Tab20, LZZ+19]. Early and
reliable precipitation forecasting is thus critical for regional water
resource management, which increasingly impacts large swaths of the population
[AFH+15]. Many studies have shown that the sea surface temperature (SST) over
various regions of the ocean, such as the El Niño-Southern Oscillation (ENSO),
are predictive of precipitation in the United States [MYR+18, Dai13, LGM20].
Understanding which ocean regions are most predictive is challenging, however,
due to a short observational record and strong correlations among SSTs
[SWM+21].

###  6.1 Importance of Ocean Climate Indices

We estimate the importance of different ocean regions for seasonal
precipitation forecasting using our lazy training method. The response is the
average winter precipitation over the Southwestern US, and as predictors we
use 10 ocean climate indices (OCIs), which are defined as the average
detrended SST anomalies over different ocean regions [CMA+16]. As data, we use
simulations from the Community Earth System Model-Large Ensemble project
(CESM-LENS; [KDP+15, dBS+19]). Details about data processing can be found in
Section B.6.

There are strong correlations among the various OCIs (Figure 4)) — in
particular, the various Niño indices appear to be nearly collinear. Because of
this, we would expect methods like linear regression to inaccurately estimate
coefficients and their importance (see appendix for more discussion).

We apply LazyVI to this problem by first training a two-layer neural network
with a hidden width of 50 and then removing each climate index and linearly
estimating the correction. When comparing with the dropout and retraining VI
estimates, we see that dropout drastically overestimates VI of Niño 3 and Niño
3.4 relative to retraining, and that LazyVI results in estimates much closer
to the retraining estimates. These results are consistent with recent
literature indicating that the predictive ability of Niño is often overstated
relative to other OCIs [MYR+18], suggesting that LazyVI could potentially help
us better understand the relative importance of different climate mechanisms.

Figure 4: Left: sample covariance matrix of the OCIs across the 40 LENS
ensemble members; Right: estimated VI for each OCI across 10 different
train/test splits. Major discrepancies between Dropout and LazyVI/Retrain
occur in regions of high correlation.

###  6.2 High-dimensional Seasonal Forecasting

Aggregating climate regions into OCIs is standard in the climate literature
and a critical tool for understanding climate dynamics. However, while more
difficult to interpret and estimate, disaggregating OCIs and investigating
individual SST locations offers important insights into the rapidly changing
climate system [SWM+21]. Neural networks have increasingly been used to make
these types of high-dimensional forecasts, and with that comes an increased
interest in explainability [MEUB22]. However, standard gradient-based
attribution/saliency methods used to interpret NNs, while powerful for
particular networks, are often subjective and difficult to interpret [MBEU22].

The ROAR (RemOve and Retrain) framework introduced by [H+19] offers a helpful
way to evaluate such importance measures. This work provides a retraining-
based benchmark for evaluating NN attribution/saliency methods by removing
variables in order of estimated importance and measuring the drop in
predictive power. This work finds that many common attribution methods are no
more informative than a random baseline, and aruges that retraining the
network after dropping out variables is key in understanding this behavior.

Using all summer SSTs across the Pacific basin on a $10^{\circ}\times
10^{\circ}$ grid, for a total of 220 predictors, we show that LazyVI can
achieve similar results to retraining in the ROAR framework at a computational
speed-up. We train fully-connected three-layer neural network of widths (100,
50) on all variables, and from this trained network, we estimate feature
importance using the baseline Gradients importance method (GRAD, [SVZ13]). We
then remove $t=(.1,.25,.5,.75,.9,.99)$ proportion of variables by removing
them in order of GRAD importance (in addition to a random ordering as a
baseline, see Figure 11 in Section B.6 for a visualization of this procedure).
We estimate model performance on these modified datasets using the Dropout,
Retraining (ROAR) and LazyVI approaches and find that LazyVI closely
approximates ROAR in nearly half the time (Figure 5) and the ordering of
variables removed does not matter much; Dropout, on the other hand, vastly
overestimates the degradation of the model performance, even with a small
number of variables removed.

Figure 5:  Average MSE across 5 runs after removing increasing proportions of
variables with the specified importance orderings (see Figure 11 for more
detail). On average, LazyVI took 1.8s and Retrain took 3.2s.

##  7 Discussion and Extended Applications

Assessing variable importance in machine learning is a vital task as learning-
based tools are increasingly integrated into societally-impactful systems,
including autonomous vehicles, financial and healthcare decision-making, and
social and criminal justice. In this work, we propose a method, LazyVI, for
efficiently estimating variable importance based on a linearization of a fully
trained neural network. We prove that our method provides an accurate estimate
of VI and can achieve the same rate of accuracy as a computationally expensive
retraining method nearly as quickly as the inaccurate dropout method. We
further show how to construct confidence intervals around these estimates.

The theory developed in this paper provides an important step toward making
interpretability in neural networks more computationally efficient, and we
suspect this theoretical framework can extend to other settings, which we
discuss here.

###  7.1 Early Stopping and Regularization

A potential alternative to our proposed LazyVI method is to first train a full
model (as we do) and then train the reduced model using a gradient-based
method initialized with the full model parameters and stopped early. Empirical
evidence suggests that this approach would have similar speed and accuracy to
our LazyVI approach due to the implicit regularization associated with early
stopping. This approach has the potential to extend LazyVI to far more
complicated architectures than the standard feedforward networks we have
experimented with thus far. As a proof of concept, we train a convolutional
neural network on the MNIST benchark dataset, and then follow the ROAR
procedure with random ordering. We see in Figure 6 that Dropout results in a
consistent decline in predictive performance when variables are removed, while
the accuracy remains relatively high when the network is retrained up until
around 75% of variables are removed - remarkable, given how relatively
uninformative the image appears (Figure 6, bottom). Importantly, we see that
taking a single step from a model initialized at the full model parameters
(LazyVI-ES) results in nearly identical performance to the full retraining at
a $5\times$ speed-up.

Figure 6:  (Left) Average accuracy of estimator across 5 repetitions. (Right)
Average computation time for a single run of each method. (Bottom) Example
MNIST image with different proportions of variables removed. Early stopping is
able to nearly identically mimic a full model retrain at a fraction of the
time.

While these results are promising, we currently lack theoretical guarantees
for early stopping in this setting. It is possible that our theoretical
results could lead to new insights into early stopping for assessing VI due to
the intimate connection between kernel ridge regression and early stopping
algorithms [RWY14]. In fact, if the eigenvalues of the NTK matrix at the full
model initialization decay in a sufficiently fast rate, early stopping of the
reduced model training should give an as good estimate of the reduced model.
However, analyzing early stopping in this setting requires characterizing the
spectrum of the NTK with _the full model initialization_ , whereas most
spectral properties of the NTK have been developed under the assumption of a
random initialization [NMM21, MZ20]. Better understanding the NTK spectrum
after full-model initialization in the future could provide new insights into
fast algorithms for VI estimation.

###  7.2 Shapley Values

When features are correlated, the quantity VI defined in (7) tends to zero.
Recent work proposes using Shapley values to measure variable importance,
arguing that their handling of correlated variables, which assigns similar
positive weights to correlated important variables, is desirable [OP16, WF20]
in some settings. These papers also note that Shapley values are prohibitively
expensive to compute, as they require fitting a new model for each of the
$2^{p}$ possible subsets of variables. However, we note that computing the
Shapley values requires many calculations of the quantity in (7); an important
avenue is investigating the use of our LazyVI framework to accelerate the
computation of Shapley values. We perform a preliminary experiment on
calculating Shapley values using our LazyVI framework, and compare it with the
retraining method used in [WF20]. We estimate the Shapley values for a sparse
high-dimensional data generated by a logistic model, and perform the
retraining/lazyVI method on a two-layer neural network. When using LazyVI
training, the computation is roughly $5$ times faster and the estimated
Shapley values are close to retraining. Moreover, when the sample size is
relatively small with respective to the dimension, we observe Lazy training
has a smaller variance of estimated Shapley values on non-important variables
than retraining method, due to the regularization proposed in our method. See
Section B.7 for our initial exploration into this line of work.

## Acknowledgements

This work was supported by AFOSR FA9550-18-1-0166, DOE DE-AC02-06CH113575, NSF
OAC-1934637, NSF DMS-1930049 and NSF DMS-2023109.

## References

  * [AFH+15] Amir AghaKouchak, David Feldman, Martin Hoerling, Travis Huxman, and Jay Lund.  Water and climate: Recognize anthropogenic drought.  Nature, 524(7566), August 2015.
  * [Bar94] Andrew R Barron.  Approximation and estimation bounds for artificial neural networks.  Machine learning, 14(1):115–133, 1994.
  * [BBM+15] Sebastian Bach, Alexander Binder, Grégoire Montavon, Frederick Klauschen, Klaus-Robert Müller, and Wojciech Samek.  On pixel-wise explanations for non-linear classifier decisions by layer-wise relevance propagation.  PloS one, 10(7):e0130140, 2015.
  * [BC18] Rina Foygel Barber and Emmanuel J. Candes.  A knockoff filter for high-dimensional selective inference, 2018.
  * [BM02] Peter L Bartlett and Shahar Mendelson.  Rademacher and gaussian complexities: Risk bounds and structural results.  Journal of Machine Learning Research, 3(Nov):463–482, 2002.
  * [CFJL17] Emmanuel Candes, Yingying Fan, Lucas Janson, and Jinchi Lv.  Panning for Gold: Model-X Knockoffs for High-dimensional Controlled Variable Selection.  arXiv:1610.02351 [math, stat], December 2017.  arXiv: 1610.02351.
  * [CMA+16] Yang Chen, Douglas C. Morton, Niels Andela, Louis Giglio, and James T. Randerson.  How much global burned area can be forecast on seasonal time scales using sea surface temperatures?  Environmental Research Letters, 11(4):045001, March 2016.
  * [COB20] Lenaic Chizat, Edouard Oyallon, and Francis Bach.  On Lazy Training in Differentiable Programming.  arXiv:1812.07956, January 2020.
  * [CRG17] Chun-Hao Chang, Ladislav Rampasek, and Anna Goldenberg.  Dropout feature ranking for deep learning models.  arXiv preprint arXiv:1712.08645, 2017.
  * [Dai13] Aiguo Dai.  The influence of the inter-decadal Pacific oscillation on US precipitation during 1923–2010.  Climate Dynamics, 41(3), August 2013.
  * [dBS+19] J. de La Beaujardière, A. Banihirwe, C. F. G. Shih, K. Paul, and J. Hamman.  Ncar cesm lens cloud-optimized subset.  UCAR/NCAR Computational and Informations Systems Lab, 2019.
  * [DS95] Kjell Doksum and Alexander Samarov.  Nonparametric estimation of global functionals and a measure of the explanatory power of covariates in regression.  The Annals of Statistics, pages 1443–1473, 1995.
  * [GMR+18] Riccardo Guidotti, Anna Monreale, Salvatore Ruggieri, Franco Turini, Fosca Giannotti, and Dino Pedreschi.  A survey of methods for explaining black box models.  ACM Comput. Surv., 51(5), August 2018.
  * [H+19] Sara Hooker et al.  A benchmark for interpretability methods in deep neural networks.  In NeurIPS, 2019.
  * [HKZ12] Daniel Hsu, Sham Kakade, and Tong Zhang.  A tail inequality for quadratic forms of subgaussian random vectors.  Electronic Communications in Probability, 17:1–6, 2012.
  * [HLY19] Wei Hu, Zhiyuan Li, and Dingli Yu.  Simple and effective regularization methods for training on noisily labeled data with generalization guarantee.  arXiv preprint arXiv:1905.11368, 2019.
  * [JGH20] Arthur Jacot, Franck Gabriel, and Clément Hongler.  Neural tangent kernel: Convergence and generalization in neural networks, 2020.
  * [KDP+15] J. E. Kay, C. Deser, A. Phillips, A. Mai, C. Hannay, G. Strand, J. M. Arblaster, S. C. Bates, G. Danabasoglu, J. Edwards, M. Holland, P. Kushner, J.-F. Lamarque, D. Lawrence, K. Lindsay, A. Middleton, E. Munoz, R. Neale, K. Oleson, L. Polvani, and M. Vertenstein.  The Community Earth System Model (CESM) Large Ensemble Project: A Community Resource for Studying Climate Change in the Presence of Internal Climate Variability.  Bulletin of the American Meteorological Society, 96(8):1333–1349, August 2015.
  * [LGM20] Nathan J. L. Lenssen, Lisa Goddard, and Simon Mason.  Seasonal Forecast Skill of ENSO Teleconnection Maps.  Weather and Forecasting, 35(6):2387–2406, December 2020.  Publisher: American Meteorological Society Section: Weather and Forecasting.
  * [LGR+18] Jing Lei, Max G’Sell, Alessandro Rinaldo, Ryan J. Tibshirani, and Larry Wasserman.  Distribution-free predictive inference for regression.  Journal of the American Statistical Association, 113(523):1094–1111, 2018.
  * [LXS+20] Jaehoon Lee, Lechao Xiao, Samuel S Schoenholz, Yasaman Bahri, Roman Novak, Jascha Sohl-Dickstein, and Jeffrey Pennington.  Wide neural networks of any depth evolve as linear models under gradient descent.  Journal of Statistical Mechanics: Theory and Experiment, 2020(12):124002, Dec 2020.
  * [LZZ+19] Chao Li, Francis Zwiers, Xuebin Zhang, Gang Chen, Jian Lu, Guilong Li, Jesse Norris, Yaheng Tan, Ying Sun, and Min Liu.  Larger Increases in More Extreme Local Precipitation Events as Climate Warms.  Geophysical Research Letters, 46(12):6885–6891, 2019.
  * [MBEU22] Antonios Mamalakis, Elizabeth A. Barnes, and Imme Ebert-Uphoff.  Investigating the fidelity of explainable artificial intelligence methods for applications of convolutional neural networks in geoscience, 2022\.
  * [MEUB22] Antonios Mamalakis, Imme Ebert-Uphoff, and Elizabeth A. Barnes.  Neural network attribution methods for problems in geoscience: A novel synthetic benchmark dataset.  Environmental Data Science, 1:e8, 2022.
  * [MRT18] Mehryar Mohri, Afshin Rostamizadeh, and Ameet Talwalkar.  Foundations of machine learning.  2018\.
  * [MYR+18] Antonios Mamalakis, Jin-Yi Yu, James T. Randerson, Amir AghaKouchak, and Efi Foufoula-Georgiou.  A new interhemispheric teleconnection increases predictability of winter precipitation in southwestern US.  Nature Communications, 9(1):2332, June 2018.
  * [MZ20] Andrea Montanari and Yiqiao Zhong.  The interpolation phase transition in neural networks: Memorization and generalization under lazy training.  arXiv preprint arXiv:2007.12826, 2020.
  * [NMM21] Quynh Nguyen, Marco Mondelli, and Guido F Montufar.  Tight bounds on the smallest eigenvalue of the neural tangent kernel for deep relu networks.  In International Conference on Machine Learning, pages 8119–8129. PMLR, 2021.
  * [OP16] Art B. Owen and Clémentine Prieur.  On Shapley value for measuring importance of dependent inputs, October 2016.
  * [RR19] Cynthia Rudin and Joanna Radin.  Why are we using black box models in ai when we don’t need to? a lesson from an explainable ai competition.  Harvard Data Science Review, 1(2), 11 2019.  https://hdsr.mitpress.mit.edu/pub/f9kuryi8.
  * [RWY14] Garvesh Raskutti, Martin J. Wainwright, and Bin Yu.  Early stopping and non-parametric regression: An optimal data-dependent stopping rule.  Journal of Machine Learning Research, 2014.
  * [SGK19] Avanti Shrikumar, Peyton Greenside, and Anshul Kundaje.  Learning Important Features Through Propagating Activation Differences.  arXiv:1704.02685 [cs], October 2019.  arXiv: 1704.02685.
  * [STK+17] Daniel Smilkov, Nikhil Thorat, Been Kim, Fernanda Viégas, and Martin Wattenberg.  Smoothgrad: removing noise by adding noise.  arXiv preprint arXiv:1706.03825, 2017.
  * [STY17] Mukund Sundararajan, Ankur Taly, and Qiqi Yan.  Axiomatic attribution for deep networks.  In International Conference on Machine Learning, pages 3319–3328. PMLR, 2017.
  * [SvdLP14] Stephanie Sapp, Mark J van der Laan, and Kimberly Page.  Targeted estimation of binary variable importance measures with interval-censored outcomes.  The international journal of biostatistics, 10(1):77–97, 2014.
  * [SVZ13] Karen Simonyan, Andrea Vedaldi, and Andrew Zisserman.  Deep inside convolutional networks: Visualising image classification models and saliency maps, 2013.
  * [SWM+21] Abby Stevens, Rebecca Willett, Antonios Mamalakis, Efi Foufoula-Georgiou, Alejandro Tejedor, James T. Randerson, Padhraic Smyth, and Stephen Wright.  Graph-Guided Regularized Regression of Pacific Ocean Climate Variables to Increase Predictive Skill of Southwestern U.S. Winter Precipitation.  Journal of Climate, 34(2):737–754, January 2021.
  * [Tab20] Hossein Tabari.  Climate change impact on flood and extreme precipitation increases with water availability.  Scientific Reports, 10(1), August 2020.
  * [WF20] Brian D. Williamson and Jean Feng.  Efficient nonparametric statistical inference on population feature importance using Shapley values.  arXiv:2006.09481 [stat], June 2020.  arXiv: 2006.09481.
  * [WGSC21] Brian D. Williamson, Peter B. Gilbert, Noah R. Simon, and Marco Carone.  A general framework for inference on algorithm-agnostic variable importance.  arXiv:2004.03683, 2021.
  * [ZJ21] Lu Zhang and Lucas Janson.  Floodgate: inference for model-free variable importance.  arXiv:2007.01283 [stat], April 2021.  arXiv: 2007.01283.

##  Appendix A Supporting Lemmas and Proofs

###  A.1 Supporting Lemma

###### Assumption A.1

  * (A1)

There exists some constant $C>0$ such that, for each sequence
$f_{1},f_{2},\dots\in\mathcal{F}$ such that
$\|f_{i}-f_{0}\|_{\mathcal{F}}\rightarrow 0$,
$|V(f_{j},P_{0})-V(f_{0},P_{0})|\leq C\|f_{j}-f_{0}\|_{\mathcal{F}}^{2}$ for
each $j$ large enough;

  * (A2)

There exists some constant $\delta>0$ such that for each sequence
$\epsilon_{1},\epsilon_{2},\dots\in\mathbb{R}$ and
$h,h_{1},h_{2},\dots\in\mathbb{R}$ satisfying that $\epsilon_{j}\rightarrow 0$
and $\|h_{j}-h\|_{\infty}\rightarrow 0$, it holds that

$\sup_{f\in\mathcal{F}:\|f-f_{0}\|_{\mathcal{F}}<\delta}|\frac{V(f,P_{0}+\epsilon_{j}h_{j})-V(f,P_{0})}{\epsilon_{j}}-\dot{V}(f,P_{0};h_{j})|\rightarrow
0;$

  * (B2)

$\int[g_{n}(z)]^{2}dP_{0}(z)=o_{P}(1)$;

###### Lemma A.2

([WGSC21]) Suppose (A1-A2, B2) regularity conditions hold. Denote $f_{n}(X)$
and $f_{n,-j}(X_{-j})$ as the estimate for $f_{0}$ and $f_{0,-j}$, Then for a
predictive skill measure $V(f,P)$ satisfying conditions (A1)-(A2), (B2) in
Appendix, as long as the estimators satisfy the following condition:

$\|f_{n}-f_{0}\|_{\mathcal{F}}=O_{p}(n^{-\frac{1}{4}}),\leavevmode\nobreak\
\leavevmode\nobreak\
\|f_{n,-j}-f_{0,-j}\|_{\mathcal{F}}=O_{p}(n^{-\frac{1}{4}}),$ (29)

for all $j\in[p]$, then we have

$\begin{split}v_{n}-v_{0}=&\frac{1}{n}\sum_{i=1}^{n}\dot{V}(f_{0},P_{0};\delta_{Z_{i}}-P_{0})+O_{p}(\frac{1}{\sqrt{n}}),\\\
v_{n,-j}-v_{0,-j}=&\frac{1}{n}\sum_{i=1}^{n}\dot{V}(f_{0,-j},P_{0,-j};\delta_{Z_{i}}-P_{0,-j})+O_{p}(\frac{1}{\sqrt{n}}),\end{split}$
(30)

where $v_{n}=V(f_{n},P_{n})$ and
$v_{n,-j}=V(f_{n,-j},P_{n,-j}),\leavevmode\nobreak\ \forall j\in[p]$.

###  A.2 Missing Proofs

####  A.2.1 Proof of Lemma 4.5

Here we present the detailed proof of Lemma 4.5, which gives the empirical
estimation error bounds for the NTK kernel ridge regression estimation. The
proof follows the proof framework provided in [HLY19].

First of all, according to kernel ridge regression, denote

  * •

$\mathbf{Y}=(Y_{1},\dots,Y_{n})^{\top}$;

  * •

$\mathbf{h}_{\theta_{f}}^{(j)}=(h_{\theta_{f}}({\mathbf{X}}_{1}^{(j)}),\dots,h_{\theta_{f}}({\mathbf{X}}_{n}^{(j)}))^{\top}$;

  * •

$\mathbf{f}_{0,-j}=(f_{0,-j}({\mathbf{X}}_{1}^{(j)}),\dots,f_{0,-j}({\mathbf{X}}_{n}^{(j)}))^{\top}$.

we have

$(\tilde{h}_{\theta_{f}+\Delta\theta_{j}}({\mathbf{X}}_{1}^{(j)}),\dots,\tilde{h}_{\theta_{f}+\Delta\theta_{j}}({\mathbf{X}}_{n}^{(j)}))^{\top}={\mathbb{K}}^{(j)}({\mathbb{K}}^{(j)}+\lambda
I_{n})^{-1}(\mathbf{Y}-\mathbf{h}_{\theta_{f}}^{(j)})+\mathbf{h}_{\theta_{f}}^{(j)}.$
(31)

Recall that
$\epsilon^{(j)}=Y-\mathbb{E}(Y|X_{-j})=Y-f_{0,-j}({\mathbf{X}}^{(j)})$, we
define its observed samples as

$\bm{\epsilon}^{(j)}=\left(Y_{1}-f_{0,-j}({\mathbf{X}}_{1}^{(j)}),\dots,Y_{n}-f_{0,-j}({\mathbf{X}}_{n}^{(j)})\right)^{\top}=\mathbf{Y}-\mathbf{f}_{0,-j};$

Recall the definition of ${\mathbf{e}}^{(j)}$, we have

${\mathbf{e}}^{(j)}=\mathbf{f}_{0,-j}-\mathbf{h}_{\theta_{f}}^{(j)}.$

Hence we have

$\begin{split}&\sqrt{n}\|\tilde{h}_{\theta_{f}+\Delta\theta_{j}}(X^{(j)})-f_{0,-j}(X^{(j)})\|_{n}\\\
=&\sqrt{\sum_{i=1}^{n}\left[\tilde{h}_{\theta_{f}+\Delta\theta_{j}}({\mathbf{X}}_{i}^{(j)})-f_{0,-j}({\mathbf{X}}_{i}^{(j)})\right]^{2}}\\\
=&\|{\mathbb{K}}^{(j)}({\mathbb{K}}^{(j)}+\lambda
I_{n})^{-1}(\mathbf{Y}-\mathbf{h}_{\theta_{f}}^{(j)})+\mathbf{h}_{\theta_{f}}^{(j)}-\mathbf{f}_{0,-j}\|\\\
=&\|{\mathbb{K}}^{(j)}({\mathbb{K}}^{(j)}+\lambda
I_{n})^{-1}(\mathbf{f}_{0,-j}+\bm{\epsilon}^{(j)}-\mathbf{h}_{\theta_{f}}^{(j)})+\mathbf{h}_{\theta_{f}}^{(j)}-\mathbf{f}_{0,-j}\|\\\
=&\|{\mathbb{K}}^{(j)}({\mathbb{K}}^{(j)}+\lambda
I_{n})^{-1}({\mathbf{e}}^{(j)}+\bm{\epsilon}^{(j)})-{\mathbf{e}}^{(j)}\|\\\
=&\|{\mathbb{K}}^{(j)}({\mathbb{K}}^{(j)}+\lambda
I_{n})^{-1}\bm{\epsilon}^{(j)}-\lambda({\mathbb{K}}^{(j)}+\lambda
I_{n})^{-1}{\mathbf{e}}^{(j)}\|\\\
\leq&\|{\mathbb{K}}^{(j)}({\mathbb{K}}^{(j)}+\lambda
I_{n})^{-1}\bm{\epsilon}^{(j)}\|+\|\lambda({\mathbb{K}}^{(j)}+\lambda
I_{n})^{-1}{\mathbf{e}}^{(j)}\|.\end{split}$ (32)

According to 4.2 and [HKZ12], we have

$P\left(\|A\bm{\epsilon}^{(j)}\|^{2}/\sigma^{2}>\text{tr}(\Sigma)+2\sqrt{\text{tr}(\Sigma^{2})t}+2\|\Sigma\|t\leavevmode\nobreak\
|\leavevmode\nobreak\ {\mathbf{X}}^{(j)}\right)\leq e^{-t},$ (33)

where $A={\mathbb{K}}^{(j)}({\mathbb{K}}^{(j)}+\lambda I_{n})^{-1}$, and
$\Sigma=A^{\top}A$. Hence we have with probability at least $1-\delta$ for any
$\delta>0$, we have

$\|{\mathbb{K}}^{(j)}({\mathbb{K}}^{(j)}+\lambda
I_{n})^{-1}\bm{\epsilon}^{(j)}\|\leq\sigma\sqrt{\text{tr}(\Sigma)+2\sqrt{\text{tr}(\Sigma^{2})\log(\frac{1}{\delta})}+2\|\Sigma\|\log(\frac{1}{\delta})}.$
(34)

Let $\lambda_{1},\dots,\lambda_{n}>0$ be the eigenvalues of
${\mathbb{K}}^{(j)}$, we then have

$\begin{split}&\text{tr}[\Sigma]=\text{tr}[A^{\top}A]=\sum_{i=1}^{n}\frac{\lambda_{i}^{2}}{(\lambda_{i}+\lambda)^{2}}\leq\sum_{i=1}^{n}\frac{\lambda_{i}^{2}}{4\lambda_{i}\cdot\lambda}=\frac{\text{tr}[{\mathbb{K}}^{(j)}]}{4\lambda},\\\
&\text{tr}[\Sigma^{2}]=\text{tr}[A^{\top}A^{2}A^{\top}]=\sum_{i=1}^{n}\frac{\lambda_{i}^{4}}{(\lambda_{i}+\lambda)^{4}}\leq\sum_{i=1}^{n}\frac{\lambda_{i}^{4}}{4^{4}\lambda(\frac{\lambda_{i}}{3})^{3}}=\frac{3^{3}\text{tr}[{\mathbb{K}}^{(j)}]}{4^{4}\lambda}\leq\frac{\text{tr}[{\mathbb{K}}^{(j)}]}{4\lambda},\\\
&\|\Sigma\|=\|{\mathbb{K}}^{(j)}({\mathbb{K}}^{(j)}+\lambda
I_{n})^{-2}{\mathbb{K}}^{(j)}\|\leq 1.\end{split}$ (35)

Hence we have with probability at least $1-\delta$,

$\begin{split}&\|{\mathbb{K}}^{(j)}({\mathbb{K}}^{(j)}+\lambda
I_{n})^{-1}\bm{\epsilon}^{(j)}\|\\\
\leq&\sigma\sqrt{\frac{\text{tr}[{\mathbb{K}}^{(j)}]}{4\lambda}+2\sqrt{\frac{\text{tr}[{\mathbb{K}}^{(j)}]}{4\lambda}}+2\log(\frac{1}{\delta})}\\\
\leq&\sigma\sqrt{\frac{\text{tr}[{\mathbb{K}}^{(j)}]}{4\lambda}+2\sqrt{\frac{\text{tr}[{\mathbb{K}}^{(j)}]}{2\lambda}}+2\log(\frac{1}{\delta})}\\\
=&\sigma\sqrt{\frac{\text{tr}[{\mathbb{K}}^{(j)}]}{4\lambda}}+\sigma\sqrt{2\log(\frac{1}{\delta})}.\end{split}$
(36)

By the fact that

$\|\lambda({\mathbb{K}}^{(j)}+\lambda
I_{n})^{-1}{\mathbf{e}}^{(j)}\|=\lambda\sqrt{({\mathbf{e}}^{(j)})^{\top}({\mathbb{K}}^{(j)}+\lambda
I_{n})^{-2}{\mathbf{e}}^{(j)}},$ (37)

we have

$\|\tilde{h}_{\theta_{f}+\Delta\theta_{j}}(X^{(j)})-f_{0,-j}(X^{(j)})\|_{n}\leq\frac{{\lambda}}{\sqrt{n}}\sqrt{({\mathbf{e}}^{(j)})^{\top}[{\mathbb{K}}^{(j)}+\lambda
I_{n}]^{-2}{\mathbf{e}}^{(j)}}+\sigma\sqrt{\frac{\text{tr}[{\mathbb{K}}^{(j)}]}{4n\lambda}}+\sigma\sqrt{\frac{2}{n}\log(\frac{1}{\delta})}.$
(38)

$\square$

####  A.2.2 Equation 39 and Its Proof

Define the Hilbert norm for a function
$f(x)=\alpha^{\top}K(x,{\mathbf{X}}^{(j)}),\leavevmode\nobreak\
\forall\alpha\in\mathbb{R}^{n}$ in the NTK kernel space is:
$\|f\|_{\mathcal{H}}=\sqrt{\alpha^{\top}{\mathbb{K}}^{(j)}\alpha}$. The
following lemma is to bound the Hilbert norm for
$\tilde{h}_{\theta_{f}+\Delta\theta_{j}}-h_{\theta_{f}}$ so that we could
bound the complexity of the function class it lies in.

###### Lemma A.3

With probability at least $1-\delta$, for any $j\in[p]$ we have

$\|\tilde{h}_{\theta_{f}+\Delta\theta_{j}}(x)-h_{\theta_{f}}(x)\|_{\mathcal{H}}\leq\sqrt{({\mathbf{e}}^{(j)})^{\top}({\mathbb{K}}^{(j)}+\lambda
I_{n})^{-1}{\mathbf{e}}^{(j)}}+\frac{\sigma}{\sqrt{\lambda}}\left(\sqrt{n}+\sqrt{2\log(\frac{1}{\delta})}\right).$
(39)

[Proof.] Recall that
$\tilde{h}_{\theta_{f}+\Delta\theta_{j}}(x)=\text{ker}_{\theta_{f}}(x,{\mathbf{X}}^{(j)})({\mathbb{K}}^{(j)}+\lambda
I_{n})^{-1}(\mathbf{Y}-\mathbf{h}_{\theta_{f}}^{(j)})+h_{\theta_{f}}(x)$.
Based on the fact that
$\mathbf{Y}-\mathbf{h}_{\theta_{f}}^{(j)}={\mathbf{e}}^{(j)}+\bm{\epsilon}^{(j)}$,
we have

$\begin{split}&\|\tilde{h}_{\theta_{f}+\Delta\theta_{j}}(x)-h_{\theta_{f}}(x)\|_{\mathcal{H}}\\\
=&\|(\mathbf{Y}-\mathbf{h}_{\theta_{f}}^{(j)})^{\top}({\mathbb{K}}^{(j)}+\lambda
I_{n})^{-1}\text{ker}_{\theta_{f}}({\mathbf{X}}^{(j)},x)\|_{\mathcal{H}}\\\
=&\sqrt{({\mathbf{e}}^{(j)}+\bm{\epsilon}^{(j)})^{\top}({\mathbb{K}}^{(j)}+\lambda
I_{n})^{-1}{\mathbb{K}}^{(j)}({\mathbb{K}}^{(j)}+\lambda
I_{n})^{-1}({\mathbf{e}}^{(j)}+\bm{\epsilon}^{(j)})}\\\
\leq&\sqrt{({\mathbf{e}}^{(j)}+\bm{\epsilon}^{(j)})^{\top}({\mathbb{K}}^{(j)}+\lambda
I_{n})^{-1}({\mathbf{e}}^{(j)}+\bm{\epsilon}^{(j)})}\\\
\leq&\sqrt{({\mathbf{e}}^{(j)})^{\top}({\mathbb{K}}^{(j)}+\lambda
I_{n})^{-1}{\mathbf{e}}^{(j)}}+\sqrt{(\bm{\epsilon}^{(j)})^{\top}({\mathbb{K}}^{(j)}+\lambda
I_{n})^{-1}\bm{\epsilon}^{(j)}}\\\
\leq&\sqrt{({\mathbf{e}}^{(j)})^{\top}({\mathbb{K}}^{(j)}+\lambda
I_{n})^{-1}{\mathbf{e}}^{(j)}}+\sqrt{\frac{(\bm{\epsilon}^{(j)})^{\top}\bm{\epsilon}^{(j)}}{\lambda}}.\end{split}$
(40)

Using the concentration inequality in [HKZ12] again, we have with probability
at least $1-\delta$, we have

$\sqrt{(\bm{\epsilon}^{(j)})^{\top}\bm{\epsilon}^{(j)}}\leq\sigma\sqrt{n+2\sqrt{n\log(\frac{1}{\delta})}+2\log(\frac{1}{\delta})}\leq\sigma\left(\sqrt{n}+\sqrt{2\log(\frac{1}{\delta})}\right).$
(41)

Hence we prove Equation 39 by combining Equation 40 and Equation 41.

####  A.2.3 Generalization Error Bound and Its Proof

In the following, we will bound the generalization error based on the above
empirical error bound.

###### Lemma A.4

For any $j\in[p]$, let $\|\cdot\|$ be the $L_{2}(P_{0})$ norm defined as
$\|f\|=\sqrt{\int|f(x^{(j)})|dP_{0}(x)}$, then we have with probability at
least $1-\delta$ for any $\delta>0$,

$\begin{split}\|\tilde{h}_{\theta_{f}+\Delta\theta_{j}}-f_{0,-j}\|\leq&\left\\{\frac{{\lambda}}{\sqrt{n}}\sqrt{({\mathbf{e}}^{(j)})^{\top}[{\mathbb{K}}^{(j)}+\lambda
I_{n}]^{-2}{\mathbf{e}}^{(j)}}+\sigma\sqrt{\frac{\text{tr}[{\mathbb{K}}^{(j)}]}{4n\lambda}}+\sigma\sqrt{\frac{2}{n}\log(\frac{3}{\delta})}\right\\}\\\
&+\frac{2\sqrt{\text{tr}[{\mathbb{K}}^{(j)}]}}{n}\left[O(1)+\frac{\sigma}{\sqrt{\lambda}}(\sqrt{n}+\sqrt{2\log(3/\delta)})\right]+\sqrt{\frac{\log(3/\delta)}{2n}}.\\\
\end{split}$ (42)

Under 4.1 and 4.2, when we take the penalty parameter in the rate
$\lambda=O(\sqrt{n})$, we have with high probability that
$\|\tilde{h}_{\theta_{f}+\Delta\theta_{j}}-f_{0,-j}\|\leq O_{p}(n^{-1/4})$.

[Proof.] According to Lemma 4.5, we know that with probability at least
$1-\delta/3$,

$\|\tilde{h}_{\theta_{f}+\Delta\theta_{j}}-f_{0,-j}\|_{n}\leq\frac{{\lambda}}{\sqrt{n}}\sqrt{({\mathbf{e}}^{(j)})^{\top}[{\mathbb{K}}^{(j)}+\lambda
I_{n}]^{-2}{\mathbf{e}}^{(j)}}+\sigma\sqrt{\frac{\text{tr}[{\mathbb{K}}^{(j)}]}{4n\lambda}}+\sigma\sqrt{\frac{2}{n}\log(\frac{3}{\delta})}.$
(43)

By [BM02], we know that the empirical Rademacher complexity for a function
class
$\mathcal{F}_{B}=\\{f(x)=\alpha^{\top}\text{ker}_{\theta_{f}}({\mathbf{X}}^{(j)},x):\|f\|_{\mathcal{H}}\leq
B\\}$ is bounded as

$\hat{\mathcal{R}}_{S}(\mathcal{F}_{B})\leq\frac{B\sqrt{\text{tr}[{\mathbb{K}}^{(j)}]}}{n}.$

According to [MRT18], with probability at least $1-\delta/3$, we have

$\begin{split}&\sup_{\tilde{h}_{\theta_{f}+\Delta\theta_{j}}-h_{\theta_{f}}\in\mathcal{F}}\left\\{\left\|\tilde{h}_{\theta_{f}+\Delta\theta_{j}}(x^{(j)})-h_{\theta_{f}}(x^{(j)})-\left(f_{0,-j}(x^{(j)})-h_{\theta_{f}}(x^{(j)})\right)\right\|-\|\tilde{h}_{\theta_{f}+\Delta\theta_{j}}(x)-f_{0,-j}(x^{(j)})\|_{n}\right\\}\\\
\leq&2\hat{\mathcal{R}}_{S}(\mathcal{F})+\sqrt{\frac{\log(3/\delta)}{2n}}.\end{split}$
(44)

From 4.1 and Equation 39, we have with probability $1-\delta/3$

$\|\tilde{h}_{\theta_{f}+\Delta\theta_{j}}(x)-h_{\theta_{f}}(x)\|_{\mathcal{H}}:=B^{\prime}\leq
O(1)+\frac{\sigma}{\sqrt{\lambda}}\left(\sqrt{n}+\sqrt{2\log(\frac{3}{\delta})}\right).$
(45)

Then we have with probability $1-\delta$,

$\begin{split}&\|\tilde{h}_{\theta_{f}+\Delta\theta_{j}}-f_{0,-j}\|\\\
\leq&\|\tilde{h}_{\theta_{f}+\Delta\theta_{j}}-f_{0,-j}\|_{n}+2\hat{\mathcal{R}}_{S}(\mathcal{F})+\sqrt{\frac{\log(3/\delta)}{2n}}\\\
\leq&\left\\{\frac{{\lambda}}{\sqrt{n}}\sqrt{({\mathbf{e}}^{(j)})^{\top}[{\mathbb{K}}^{(j)}+\lambda
I_{n}]^{-2}{\mathbf{e}}^{(j)}}+\sigma\sqrt{\frac{\text{tr}[{\mathbb{K}}^{(j)}]}{4n\lambda}}+\sigma\sqrt{\frac{2}{n}\log(\frac{3}{\delta})}\right\\}+\frac{2B^{\prime}\sqrt{\text{tr}[{\mathbb{K}}^{(j)}]}}{n}+\sqrt{\frac{\log(3/\delta)}{2n}}\\\
\leq&\left\\{\frac{{\lambda}}{\sqrt{n}}\sqrt{({\mathbf{e}}^{(j)})^{\top}[{\mathbb{K}}^{(j)}+\lambda
I_{n}]^{-2}{\mathbf{e}}^{(j)}}+\sigma\sqrt{\frac{\text{tr}[{\mathbb{K}}^{(j)}]}{4n\lambda}}+\sigma\sqrt{\frac{2}{n}\log(\frac{3}{\delta})}\right\\}\\\
&+\frac{2\sqrt{\text{tr}[{\mathbb{K}}^{(j)}]}}{n}\left[O_{p}(1)+\frac{\sigma}{\sqrt{\lambda}}(\sqrt{n}+\sqrt{2\log(3/\delta)})\right]+\sqrt{\frac{\log(3/\delta)}{2n}}\\\
\end{split}$ (46)

By the assumptions that $\|[{\mathbb{K}}^{(j)}+\lambda
I_{n}]^{-1}{\mathbf{e}}^{(j)}\|^{2}=O_{p}(1/\sqrt{n})$ and
$\text{tr}[{\mathbb{K}}^{(j)}]=O_{p}(n)$ in 4.1 (a) (b), when we take
$\lambda=O(\sqrt{n})$, we have

$\|\tilde{h}_{\theta_{f}+\Delta\theta_{j}}-f_{0,-j}\|\leq O_{p}(n^{-1/4}).$
(47)

####  A.2.4 Proof of Lemma 4.6

Lemma 4.6 For a large neural network whose width is in the order of
$O(\sqrt{n})$ where $n$ is the training sample size, our lazy trained neural
network is close to its linearization with high probability:

$\|\tilde{h}_{\theta_{f}+\Delta\theta_{j}}(x)-h_{\theta_{f}+\Delta\theta_{j}}(x)\|\leq
O(n^{-1/4}).$ (48)

[Proof.] Since
$\tilde{h}_{\theta_{f}+\Delta\theta_{j}}(x)=h_{\theta_{f}}+\Delta\theta_{j}^{\top}\nabla_{\theta}h_{\theta}(x)|_{\theta=\theta_{f}}$
is a linearization of $h_{\theta_{f}+\Delta\theta_{j}}(x)$ around the
initialization $\theta_{f}$, according to Theorem 2.1 in [LXS+20], when the
neural network has a width $M$, the neural network is close to its
linearization with probability arbitrarily close to 1:

$\|\tilde{h}_{\theta_{f}+\Delta\theta_{j}}(x)-h_{\theta_{f}+\Delta\theta_{j}}(x)\|_{2}=O(\frac{1}{\sqrt{M}}).$
(49)

Specifically, when the neural network $M$ takes the order of $O(\sqrt{n})$, we
have
$\|\tilde{h}_{\theta_{f}+\Delta\theta_{j}}(x)-h_{\theta_{f}+\Delta\theta_{j}}(x)\|_{2}=O(n^{-1/4}).$

####  A.2.5 Proof of the Main Theorem (Theorem 4.4)

Based on Lemma A.4 and Lemma 4.6 , for a neural network with width at least
$M=O(\sqrt{n})$ when the assumptions hold true, by triangular inequality we
have

$\|h_{\theta_{f}+\Delta\theta_{j}}-f_{0,-j}\|\leq\|h_{\theta_{f}+\Delta\theta_{j}}-\tilde{h}_{\theta_{f}+\Delta\theta_{j}}\|+\|\tilde{h}_{\theta_{f}+\Delta\theta_{j}}-f_{0,-j}\|=O_{p}(n^{-1/4}).$
(50)

This holds true for any $j\in[p]$. Then by Lemma A.2, we finish the proof for
Theorem 4.4.

####  A.2.6 Proof of 3.1

The density of $X_{1}$ given $X_{2}$ in the setting of 3.1 is:

$\begin{split}f(x_{1}|x_{2};\rho,\sigma)=\frac{1}{2\pi\sigma^{2}\sqrt{1-\rho^{2}}}\exp{\left\\{-\frac{1}{2(1-\rho^{2})\sigma^{2}}(x_{1}^{2}-2\rho
x_{1}x_{2}+x_{2}^{2})\right\\}},\end{split}$ (51)

thus we have $X_{1}|X_{2}\sim\mathcal{N}(\rho X_{2},(1-\rho^{2})\sigma^{2})$.

##  Appendix B Additional Experiments

###  B.1 Trace Divergence Rate of the Neural Tangent Kernel Matrix

In 4.1(b), we assume the trace of the neural tangent kernel matrix with full-
model parameters as initialization diverges in the order of $n$ in
probability: $\text{tr}({\mathbb{K}}^{(j)})=O_{p}(n)$. In the following
experiment, we’ll verify this through a simulation.

We consider a two-layer neural network with 128 nodes in the hidden layer. The
data is generated from a sparse linear model
$Y=1.5X_{1}+1.2X_{2}+X_{3}+\epsilon$ where
$\epsilon\sim\mathcal{N}(0,0.1^{2})$ and we have $6$ predictors
$X_{1},X_{2},\dots,X_{6}$ generated from a normal distribution
$\mathcal{N}(0,I_{6}+C)$ with $C_{1,2}=C_{2,1}=0.5$ and $C_{i,j}=0$ for
$(i,j)\notin\\{(1,2),(2,1)\\}$ (this is to add some correlation to the
predictors); The total sample size of the data varies in the set
$\\{1000,1100,\dots,4000\\}$; Among these simulated samples at each sample
size, $2/3$ data are sampled into the training set and $1/3$ samples fall in
the testing set.

We’ll first train the full NN model on the training set and get inferred
parameters $\hat{\theta}_{f}$ in the neural network from the full training
data. Then we use $\hat{\theta}_{f}$ as initialization for the reduced neural
network. Then we could calculate the neural tangent kernel matrix and its
corresponding trace the testing data with one feature (e.g. the first feature)
dropped and replaced with $0$ (its population mean). We will repeat this
process at each sample size level five times (with different random seed to
generate data each time), and record the traces with respect to the test data
sample size. As shown in Figure 7, there is a clear trend that the trace
diverges linearly as the sample size, which numerically verifies of 4.1 (b).

Figure 7: The divergence of the NTK matrix trace with respect to the sample
size
We calculate the Neural Tangent Kernel(NTK) matrix from a two-layer neural
network with $128$ hidden nodes. The initialization of the network is the full
model parameters trained by the training data, and the input is the test data
with the first feature replaced with $0$s. The trace of this NTK matrix is
plotted against the sample size of the test data. There is a clear linear
pattern in the divergence of the NTK matrix trace w.r.t. the sample size,
which verifies Assumption (b).

###  B.2 Choosing the Regularization Parameter

LazyVI involves solving a ridge regression to estimate the difference between
the full and reduced model parameters. For variable $j$, we choose the
regularization parameter $\lambda_{j}$ through K-fold cross validation on the
prediction made using the estimated $\Delta\theta_{j}^{\lambda}$. Algorithm 2
below shows the entire procedure.

Algorithm 2 K-Fold CV for $\lambda_{j}$

0:
$\\{{\mathbf{X}}_{i}^{(j)},Y_{i},{\mathbf{e}}_{i}^{(j)},\Phi_{i}^{(j)}\\}_{i=1}^{n_{1}}$
and $\theta_{f}$ from Algorithm 1 in main paper; candidate $\lambda$ values
$\Lambda$

Partition $[n_{1}]$ into $K$ subsets, each denoted $S_{t}$

for $\lambda\in\Lambda$ do

for $k=1,\dots,K$ do

$\Delta\theta_{j}^{\lambda}=\arg\min_{\omega\in\mathbb{R}^{M}}\frac{1}{n_{1}-|S_{k}|}\sum_{i\notin
S_{k}}[{\mathbf{e}}_{i}^{(j)}-\langle\omega,\Phi_{i}^{(j)}\rangle]^{2}+\lambda\|\omega\|_{2}^{2}$

$\hat{Y}_{i}=h_{\theta_{f}+\Delta\theta_{j}^{\lambda}}({\mathbf{X}}_{i}^{(j)})$
for $i\in S_{k}$

$\epsilon_{\lambda,k}=\frac{1}{|S_{k}|}\sum_{i\in
S_{k}}(Y_{i}-\hat{Y}_{i})^{2}$

end for

$\epsilon_{\lambda}=\frac{1}{K}\sum_{k=1}^{K}\epsilon_{\lambda,k}$

end for

$\hat{\lambda}_{j}=\mathop{\rm
arg\min}_{\lambda}\\{\epsilon_{\lambda}\\}_{\lambda\in\Lambda}$

###  B.3 Full Linear Experiment

Figure 8 shows the distribution of computation time vs. VI estimation accuracy
for three different groups of variables (important and correlated, important
and uncorrelated, unimportant and uncorrelated). We see that LazyVI and
retrain are both accurate across all groups of variables, but LazyVI is much
faster. Dropout is consistently the fastest method, but is highly inaccurate
in estimating vi for the first group of variables due to their strong
correlations. Also in Figure 8 we show the empirical coverage of the LazyVI
and retrain 95% confidence intervals. We see that both retrain and LazyVI
achieve desirable coverage for the three important variables; poor coverage of
unimportant variables is expected and possibly remedied with a sample-
splitting procedure [WGSC21].

Figure 8: Distribution of computation time vs. estimation error relative to
retrain ($\hat{\text{\sc vi}}-\widehat{\text{\sc vi}}^{(\text{\sc rt})})$ for
three different groups of variables: important, correlated
$(\\{X_{1},X_{2}\\})$; important, uncorrelated ($X_{3}$); and unimportant,
uncorrelated $(\\{X_{4},X_{5},X_{6}\\})$. 2D box plots show quantiles across
10 repetitions.

###  B.4 Impact of Lazy Initialization

As discussed in the main paper, the initialization of the LazyVI procedure
plays a significant role in the accuracy of its estimates. Figure 9 shows the
distribution of the VI error for dropout, LazyVI with a good initialization,
and LazyVI with a random initialization across 10 repetitions. We see that the
random initialization results in less accurate estimates with high variance.

Figure 9:  Distribution of $\text{VI}-\hat{\text{VI}}$ for the first 3
variables for dropout, LazyVI initialized with the parameters from the full
model, and LazyVI with a random initialization.

###  B.5 Width of Training Network

Theorem 4.4 implies that LazyVI will perform well when the training network is
sufficiently wide. Figure 10 shows the empirical coverage of the 95%
confidence intervals defined in (25) (across 40 repetitions) for increasing
hidden layer widths. We see that coverage increases as the width of the
network increases, but the trade-off is that the computation time for LazyVI
also increases with the network width (although remains much faster than
retraining).

Figure 10: Left: empirical coverage of 95% confidence intervals (across 50
repetitions) of LazyVI estimates for increasing widths of the training network
for the important variables. Dotted line shows 95% coverage; Right: average
computation time for LazyVI with increasing network widths.

###  B.6 Additional Details for Seasonal Forecasting Experiment

For our real data seasonal precipitation forecasting experiment, we use
simulations from the Community Earth System Model-Large Ensemble project
(CESM-LENS; [KDP+15, dBS+19]). CESM-LENS is a 40-member ensemble of climate
simulations, where the ensemble members all have the same physics but
different initial conditions. From this dataset, we extracted monthly sea
surface temperature (SST) records from 1940-2005 on a $1.25^{\circ}\times
0.9^{\circ}$ grid. We compute SST anomalies at each grid point relative to the
time period 1950-1989444https://climatedataguide.ucar.edu/climate-data/nino-
sst-indices-nino-12-3-34-4-oni-and-tni by subtracting the monthly mean and
dividing by the monthly standard deviation, and then we linearly detrend each
time series.

To compute the 10 ocean climate indices (OCI) used in our experiment, we find
the average summer (July-October) monthly SST values of these detrended SST
anomalies over specified ocean regions. These regions are well established in
the literature; we refer to the supplement from [CMA+16] to define the
boundaries of all OCIs besides NZI, for which we use [MYR+18]. See 1 for the
specific boundaries. As a response, we use the average winter (November-March)
precipitation over part of the southwestern US (see [SWM+21]). We are
interested in predicting winter precipitation from the previous summer’s SSTs.

Ocean | OCI | Latitude | Longitude
---|---|---|---
Pacific | Niño1+2 |  $10^{\circ}$S - $0^{\circ}$ |  $90^{\circ}$W - $80^{\circ}$W
| Niño3 |  $5^{\circ}$S - $5^{\circ}$N |  $150^{\circ}$W - $90^{\circ}$W
| Niño3.4 |  $5^{\circ}$S - $5^{\circ}$N |  $170^{\circ}$W - $120^{\circ}$W
| Niño4 |  $5^{\circ}$S - $5^{\circ}$N |  $160^{\circ}$E - $150^{\circ}$W
| NZI |  $40^{\circ}$S - $25^{\circ}$S |  $170^{\circ}$E - $160^{\circ}$W
Atlantic | TNA |  $5^{\circ}$N - $25^{\circ}$N |  $55^{\circ}$W - $15^{\circ}$W
| TSA |  $20^{\circ}$S- $0^{\circ}$ |  $30^{\circ}$W - $10^{\circ}$E
Indian | SWIO |  $32^{\circ}$S - $25^{\circ}$S |  $31^{\circ}$E - $45^{\circ}$E
| WTIO |  $10^{\circ}$S- $10^{\circ}$N |  $50^{\circ}$E - $70^{\circ}$E
| SETIO |  $10^{\circ}$S- $0^{\circ}$ |  $90^{\circ}$E - $110^{\circ}$E
Table 1: Ocean climate indices (OCIs) are defined as the average of the
detrended SST anomalies across the regions indicated above.

Below, we provide an example of SST removal for the experiments in Section
6.2.

Figure 11:  For an example year (1940), Figure 11 shows the original image of
SST anomalies, and then the modified images after dropping out 50% and 75% of
the data according to the GRAD importance measure and a baseline random order.
In Figure 5, we report the decline in predictive performance when removing
variables according to these orderings.

###  B.7 Shapley Value Calculation using Lazy Training Method

As we have discussed in the paper, our method may also provide a faster
alternative when calculating Shapley values using large neural networks. We’ll
indicate this in the following experiment.

We define population Shapley values in alignment with [WF20] with an arbitrary
measure of predictiveness $V$. For any variable $X_{j}$, its Shapley value
$\psi_{j}$ is defined as the average gain in oracle predictiveness from
including feature $X_{j}$ over all possible subsets:

$\psi_{j}:=\sum_{s\in[p]\setminus\\{j\\}}\frac{1}{p}{{p-1}\choose|s|}^{-1}\left[V(f_{0,s\cup\\{j\\}},P_{0})-V(f_{0,s},P_{0})\right],$
(52)

where $P_{0}$ is the true distribution and $f_{0,s},f_{0,s\cup\\{j\\}}$ are
the oracle prediction functions over the function subset
$\mathcal{F}_{s}:=\\{f\in\mathcal{F}:f(u)=f(v)\text{ for all
}u,v\in\mathcal{X}\text{ and }u_{s}=v_{s}\\}$. Shapley values have nice
properties such as non-negativity, additivity, symmetry, zero for null
features, etc. Moreover, such defined Shapley values can assign positive
values to collinear variables that are each marginally predictive, whereas
previously defined population VIs would assign zero to all collinear
variables.

We use the same subset sampling scheme as [WF20]. However, instead of
retraining the model for each subset of features, we generalize our proposed
LazyVI method to make the computation more scalable and much faster.

We are dealing with a Logistic Model with high dimensional sparse features.
Specifically, we have $100$ features from a $\mathcal{N}(0,\Sigma_{100\times
100})$, where the variables are independent except
$\text{Corr}(X_{1},X_{2})=0.75$. The responses are binary, generated from a
logistic model: $\log\frac{\mathbb{P}(Y=1)}{1-\mathbb{P}(Y=1)}=X\beta$, where
$\beta=(5,4,3,2,1,0,\dots,0)^{\top}\in\mathbb{R}^{100}$.

We use a two-layer neural network (with $128$ hidden nodes) and the subset
sampling scheme proposed by [WF20] when calculating Shapley values. We compare
the estimated Shapley values and the computing times when we use the
retraining method (as is used in [WF20]) and the lazy training method we
proposed. We calculate the Shapley values on $20$ simulated datasets, each
with a sample size $750$. We split each dataset into a training set and test
set with sample sizes $500$ and $250$ respectively. We train all the models on
the training set, then evaluate the predictiveness loss and calculate the
Shapley values on the test set. We set the penalty parameter as $50$ in
accordance to the assumption $\lambda=O(\sqrt{n})$.

In the retraining method we’ll reconstruct two-layer neural networks with 128
hidden nodes for each subset of features; in the lazy training method however,
after training a two-layer neural network with all the features, we use
Algorithm 1 to train the new model on each subset of features (with all
features not included in the subset set equal to their mean).

| Retrain | Lazy
---|---|---
Time | 272.75s | 50.27s
Std. | 10.8s | 1.6s
Table 2: Average time to compute the Shapley values in one data

Table 2 gives the average time to calculate the Shapley values for one data
set. Lazy training speeds up the calculation by more than $5$ times. In the
meanwhile, we don’t sacrifice too much on the Shapley value estimation
performance. As shown in Figure 12, the Shapley value calculated by lazy
training is generally close to the retraining results; for those unimportant
variables (Feature ID $\geq 6$), Shapley values estimated from the retraining
method have a larger variance, as the sample size is relatively small
(training size $n=500$ while $p=100$), due to the benefit of the
regularization step we have in the LazyVI method.

Figure 12: Shapley Values calculated by Retraining vs LazyVI using a two-layer
neural network. The data is generated from a logistic model, where only the
top 5 features have non-zero weights and all the remaining features have zero
weights (The features with zero weights have Shapley values of zero). The
experiment is repeated 20 times. The colored bars are the averaged estimated
Shapley value of each feature using different methods, and the gray lines
indicate the standard deviations. We can see that for non-zero Shapley value
variables, the LazyVI estimation is close to Retraining estimation, while for
variables with zero Shapley values, retraining estimation has a larger
variance than LazyVI.