Below we address your concerns one by one:

Experimental Designs Or Analyses: Our experiments are conducted on ImageNet (3224224). We also provide the detail optimization algorithm in the third reviewer's rebuttal section and the score network training algorithm at the end of this rebuttal.

Questions For Authors:

1. While training a diffusion model for an entire dataset can be expensive, our approach offers flexibility. We can train a diffusion model on a single image (3224224), which yields a score network that is sufficient for calculating pixel-wise mutual information for that image. In our experiments, generating an attribution map for one image takes about 4 minutes on a single GTX1080Ti. Moreover, when a pretrained diffusion model is available, the computation time is greatly reduced. Once the diffusion model is available (whether pretrained or efficiently trained on a single image), the attribution procedure itself mainly involves evaluating the score function and performing numerical integration (e.g., for mutual information). Our method employs closed-form expressions for the Variance Exploding (VE) and Variance Preserving (VP) processes, allowing efficient computation. These operations scale linearly with the input size and are highly parallelizable on modern hardware.

2. Our approach does have limitations. If a pretrained diffusion model is not available, one must train a model on a per-image basis—which, while less demanding than training on a full dataset, still requires significant computational resources (around 4 minutes per image on a GTX1080Ti). Although our method generally outperforms traditional attribution methods such as DeepLIFT, InputIBA, and IBA on metrics (as shown in our experiments and in the second reviewer's rebuttal), there are scenarios where the additional computational overhead may not justify the performance gains compared to simpler methods.

3. We acknowledge that our approach may be sensitive to the specific SDE form and noise schedule used. In the revised manuscript, we will include a sensitivity analysis to show how different choices of SDE parameters (such as varying $\sigma$ in Equation (36)) affect the resulting attributions. We will also discuss potential failure cases and limitations of our framework to provide a balanced perspective. We appreciate the suggestion to extend Figure 2 with additional examples. We will supplement our current experimental section with further qualitative and quantitative evaluations, including more extensive visual comparisons and analyses across diverse datasets and network architectures.

Algorithm: Training Score Network for a Single Image

Input:

• $\boldsymbol{x}_0$: Original clean image
• $s_\theta$: Score network with parameters $\theta$
• $T$: Maximum diffusion time
• $N$: Number of training iterations

Function TrainScoreNetwork:

1. Initialize noise schedule:

   • Define $\sigma_{min}$ and $\sigma_{max}$ (e.g., 0.01 and 50)
   • Define noise schedule function $\sigma(t) = \sigma_{min} \cdot (\sigma_{max}/\sigma_{min})^t$ for $t \in [0,1]$

2. For $i = 1$ to $N$:

   • Sample random time $t \sim \mathcal{U}(0, 1)$

   • Compute noise level $\sigma_t = \sigma(t)$

   • Sample noise $\boldsymbol{\epsilon} \sim \mathcal{N}(0, \boldsymbol{I})$

   • Create noisy image $\boldsymbol{x}_t = \boldsymbol{x}_0 + \sigma_t \cdot \boldsymbol{\epsilon}$

   • Compute true score $\nabla_{\boldsymbol{x}_t} \log p(\boldsymbol{x}_t|\boldsymbol{x}_0) = -\boldsymbol{\epsilon}/\sigma_t$

   • Predicted score $\hat{\boldsymbol{s}}_\theta(\boldsymbol{x}_t, t) = s_\theta(\boldsymbol{x}_t, t)$

   • Compute loss $\mathcal{L}(\theta) = \sigma_t^2 \cdot \|\hat{\boldsymbol{s}}_\theta(\boldsymbol{x}_t, t) - (-\boldsymbol{\epsilon}/\sigma_t)\|_2^2$

   • Update $\theta$ using $\nabla_\theta \mathcal{L}(\theta)$

3. Return $s_\theta$

Below, we address your specific questions and concerns:

Theoretical Claims: Please see in 'Claims And Evidence' part of in the first reviewer’s rebuttal section.

Experimental Designs Or Analyses: Our experiments are conducted on ImageNet, as it is a standard dataset frequently employed in prior work on attribution (e.g., IBA[Schulz et al., 2020]). For the classifier, we use VGG16, which is widely used in attribution studies [Schulz et al., 2020, Zhang et al., 2021]

Questions For Authors:

1. How is the loss in Equation (33) utilized?
   In our framework, the loss in Equation (33) is designed to approximate the mutual information between the perturbed input $X_t$ and the original input $X_0$. Specifically, since increasing the noise level $t$ reduces the mutual information, we use $-t$ as a proxy for the mutual information loss ($L_{\mathrm{MI}}$). This term is then combined with the standard cross-entropy loss ($L_{\mathrm{CE}}$) for classification, yielding a total loss of the form

   $$L = L_{\mathrm{MI}} + L_{\mathrm{CE}}.$$

   Minimizing this combined loss encourages the model to find perturbations that significantly reduce mutual information—thus filtering out uninformative features—while ensuring that the classifier's predictions remain correct.

2. How is mutual information used to perform feature attribution?
   Our method leverages the dynamics of mutual information along a continuous-time perturbation path defined by an SDE. By linking the time derivative of the mutual information (via its connection to the Fisher divergence) to the evolution of noise in the input, we obtain a quantitative measure of how much each feature contributes to preserving the model's output. Concretely, after complete optimization, we derive a tensor $t$ (with the same dimensions as the input) where each element indicates the “time-to-noise” required to significantly reduce the mutual information of that feature. Features with lower $t$ values are deemed more important because they are less robust to noise injection—i.e., perturbing these features causes larger changes in the mutual information and ultimately in the model’s output. Moreover, using Equations (49) and (50), we can further refine these scores by integrating the mutual information dynamics over time, thus arriving at a more precise attribution map.

We provide the detail algorithm for the optimaztion below for clarification.

Algorithm: Training $\boldsymbol{\tau}_\theta$ for Diffusion-based Feature Attribution

Input:

• $\boldsymbol{\tau}_\theta$: U-Net model with parameters $\theta$
• $f$: Pre-trained classifier
• dataset: Training dataset of images and labels
• diffusion_type: Either "VE" or "VP"

Function TrainDiffusionAttribution($\boldsymbol{\tau}_\theta$, $f$, dataset, diffusion_type):

1. Initialize diffusion parameters:

   • If diffusion_type == "VE":
     • $\sigma_{min} = 0.01$
     • $\sigma_{max} =$ appropriate value for dataset
   • If diffusion_type == "VP":
     • $\beta_{min} = 0.1$
     • $\beta_{max} = 20$

2. While not converged:

   • Sample batch $(\boldsymbol{x}, \boldsymbol{y})$ from dataset

   • $\boldsymbol{t} = \boldsymbol{\tau}_\theta(\boldsymbol{x})$

   • Sample noise $\boldsymbol{\epsilon} \sim \mathcal{N}(0, \boldsymbol{I})$

   • Generate perturbed image $\boldsymbol{z}$:

     • If diffusion_type == "VE":
       • $\sigma(\boldsymbol{t}) = \sigma_{min} \cdot (\sigma_{max}/\sigma_{min})^{\boldsymbol{t}}$
       • $\boldsymbol{z} = \boldsymbol{x} + \sigma(\boldsymbol{t}) \odot \boldsymbol{\epsilon}$
     • If diffusion_type == "VP":
       • $\beta(\boldsymbol{t}) = \beta_{min} + \boldsymbol{t} \cdot (\beta_{max} - \beta_{min})$
       • $\int_0^{\boldsymbol{t}} \beta(s) ds = \beta_{min} \cdot \boldsymbol{t} + \frac{1}{2} \cdot (\beta_{max} - \beta_{min}) \cdot \boldsymbol{t}^2$
       • $\boldsymbol{z} = \boldsymbol{x} \odot e^{-\frac{1}{2} \int_0^{\boldsymbol{t}} \beta(s) ds} + \sqrt{1 - e^{-\int_0^{\boldsymbol{t}} \beta(s) ds}} \odot \boldsymbol{\epsilon}$

   • $\hat{\boldsymbol{y}} = f(\boldsymbol{z})$

   • $\mathcal{L}_{CE} = -\sum \boldsymbol{y} \log \hat{\boldsymbol{y}}$

   • $\mathcal{L}_{MI} = \text{mean}(\boldsymbol{t})$

   • $\mathcal{L}_{MI} = \text{mean}(\boldsymbol{t})$

   • $\mathcal{L}$ = $\mathcal{L}_{CE} + \mathcal{L}_{MI}$

   • Update $\theta$ using $\nabla_\theta \mathcal{L}$

3. Return $\boldsymbol{\tau}_\theta$

Below are our responses to the specific concerns:

Other Strengths And Weaknesses:

1. Technical Contribution Clarification: While it is true that the mutual information dynamics under perturbations build on existing ideas, our work significantly extends these concepts by analyzing general SDEs (i.e., with time-dependent drift $\mu(t)$ and diffusion $\sigma(t)$). This extension is nontrivial and enables us to integrate the KL–Fisher divergence relationship into a new optimization framework for feature attribution—offering both theoretical insights and practical advantages in explaining neural network predictions.

2. Restating and Verifying Assumptions for Theorem 3.1 and Lemma 3.2: We assume that for every $t \geq 0$, the densities $p_t(\mathbf{x})$ and $q_t(\mathbf{x})$ are twice continuously differentiable and decay sufficiently fast at infinity (i.e., there exists some $m>0$ such that $\lim_{\|\mathbf{x}\|\to\infty} \|\mathbf{x}\|^m\,p_t(\mathbf{x}) = 0$, and similarly for $q_t$). These conditions ensure that all integrals are finite, boundary terms vanish under integration by parts, and the required interchanges of differentiation and integration are justified.

3. Mutual Information Loss Proxy ($L_{\mathrm{MI}}$): Recognizing that increasing $t$ (i.e., injecting more noise) monotonically reduces mutual information, we adopt $-t$ as an effective proxy for $L_{\mathrm{MI}}$. This surrogate is both computationally efficient and consistent with prior works (e.g., IBA and InputIBA), even though the theoretical derivation of the mutual information dynamics is not directly reflected in the loss formulation.

Questions For Authors:

1. Discretization Error in Mutual Information Integration:
   Once the optimization yields the tensor $t$, the integration limits for the mutual information calculation are determined, and we can choose sufficiently small step sizes for the numerical integration. This ensures that discretization errors are negligible and do not affect either the computational efficiency or the accuracy of the mutual information integration.

2. Applicability to More Complex SDEs:
   Although our current experiments use SDEs with explicit solutions (e.g., the Variance Preserving or Variance Exploding processes), our framework is general. For more complex SDEs without explicit solutions, one can rely on established numerical SDE integration methods, albeit with higher computational cost, without compromising the overall approach.

3. Comparison with InputIBA: Our experiments—particularly the Effective Heat Ratios (EHR) evaluation—demonstrate that our method clearly outperforms InputIBA. We observe a significantly higher EHR score for our approach, which indicates that our attributions are more concentrated in semantically meaningful regions. Furthermore, we provide more experimental results after the 'Other Comments Or Suggestions' sector.

4. The Range of $\sigma$ in Equation (36):
   Our approach is designed to explore the application of state-of-the-art diffusion models to feature attribution. At present, we employ the Variance Exploding (VE) and Variance Preserving (VP) processes, as these are the most representative and well-established diffusion models in the literature. In future work, we intend to investigate broader cases, including scenarios with larger ranges of $\sigma$ in Equation (36), to further generalize and enhance our method.

Other Comments Or Suggestions: We appreciate the reviewer’s note regarding the formatting issues. These will be corrected in the final version.

More Experiments

1. Bounding Boxes Evaluation [Schulz et al., 2020]

We leverage human-annotated bounding boxes from the ImageNet dataset to evaluate localization.

Box-Ratio Results:

2. Segmentation-based Ratio Evaluation Using semantic segmentation masks from the FSS-1000 dataset, we replace bounding boxes with detailed segmentation regions and compute the Segmentation-based Ratio (SR) in the same manner as the Box-Ratio. | Method | SR | |----------------------|-------------------| | IBA | 0.488 ± 0.004 | | InputIBA | 0.468 ± 0.003 | | Ours | 0.501 ± 0.003 | | Integrated Gradients | 0.079 ± 0.006 | | Guided-BP | 0.078 ± 0.009 | | Deep-Lift | 0.080 ± 0.002 | | HSIC-Attribution | 0.377 ± 0.005 |

Below, we address each of your concerns in detail.

Claims And Evidence:

Comparison with Theorem 1 in Lyu et al. (2012):
Theorem 1 in Lyu et al. (2012) considers the simple case of the SDE
$\mathrm{d}Y_t = \mathrm{d}W_t,$ where the drift is $\mu(t)=0$ and the diffusion coefficient is $\sigma(t)=1$. This setup corresponds to a basic Wiener process.

In contrast, our work extends this result to the much more general case of
$\mathrm{d}X_t = \mu(t)\ \mathrm{d}t + \sigma(t)\ \mathrm{d}W_t,$ where both the drift $\mu(t)$ and the diffusion $\sigma(t)$ are allowed to be time-dependent. This generalization is significant because:

1. Broader SDE Models:
   While Lyu et al. (2012) deal with a noise model where $\mu(t)=0$ and $\sigma(t)=1$, our formulation covers a wide range of stochastic processes. This is essential for practical applications, such as those in modern diffusion models DDPM likeand score-based diffusion models, where the noise characteristics evolve over time.
2. Handling of Additional Complexity:
   Extending the proof in Lyu et al. (2012) to accommodate non-zero drift and time-varying diffusion requires new insights. For example, our proof first builds the intuition that the drift term $\mu(t)$ can be canceled out under certain conditions—a step that is unnecessary in the Wiener process setting. Similarly, the treatment of $\sigma(t)$ involves more delicate handling due to its time dependency.

We follow the proof structure of Lyu et al. (2012) to enhance readability and make it easier for readers to identify which parts of our derivation are directly inherited from prior work and which parts constitute our innovations. We appreciate the reviewer’s suggestion to explicitly mention that our proof builds on the ideas of Lyu et al. (2012), and we will clarify this point in the revised manuscript.

Theoretical Claims:

1. Please see 'Claims And Evidence' part.
2. • The boundary terms vanish due to suitable decay at infinity or zero boundary values for $p_t$.
   • $\int p_t\ \nabla(\log p_t)\ dx = \int \nabla p_t\,dx = 0.$

Hence, $\int \nabla p_t \log p_t dx \=\ 0.$

3. We assume that for every $t\geq 0$, the densities $p_t(\mathbf{x})$ and $q_t(\mathbf{x})$ are twice continuously differentiable and decay sufficiently fast at infinity (i.e., $\lim_{\|\mathbf{x}\|\to\infty}\|\mathbf{x}\|^m\,p_t(\mathbf{x})=0$ and similarly for $q_t$ so that all integrals are finite and boundary terms vanish.
4. $\Delta$ is the Laplacian operator.

Experimental Designs Or Analyses:

1. We treat $\beta$ as a hyperparameter. In our experiments, we typically set $\beta = 1$, following the convention in prior works such as IBA and InputIBA. And we will clarify these details in the revised manuscript.
2. $L_{\mathrm{MI}}$ in our formulation leverages the observation that increasing $t$ (i.e., injecting more noise) monotonically decreases the mutual information between the original input and the perturbed input; thus, we can use $-t$ as an effective proxy for $L_{\mathrm{MI}}$.$(L_{\mathrm{CE}}$ is the standard cross-entropy loss (please see before the Equation (33)), ensuring that the classifier’s predictions remain correct for the perturbed input. We will provided the training algorithm in the following.
3. We can either train a diffusion model from scratch or use a pretrained one for efficiency. If a pretrained model is used, it must be sufficiently close to the domain of interest so that its learned score function remains accurate. it is possible to train a diffusion model on a single image, and we provide the optimization training algorithm in the third reviewer's rebuttal section and score network training algorithm in the fourth reviewer’s rebuttal section due to space constraints.
4. After the complete optimization, we obtain a tensor $t$ with the same dimensions as the input image; this tensor effectively encodes the "time-to-noise" for each pixel, which we can interpret as a preliminary measure of feature importance—pixels with lower $t$ values indicate features that are more critical for preserving the model’s output. Moreover, by further computing the integrated mutual information along the diffusion path using Equations (49) and (50), we can yield a more precise attribution map.
5. Currently, our experiments are conducted on ImageNet, as it is a standard dataset frequently employed in prior work on attribution (e.g., IBA[Schulz et al., 2020]). For the classifier, we use VGG16, which is widely used in attribution studies [Schulz et al., 2020, Zhang et al., 2021] and provides a reliable benchmark. Moreover, incorporating additional classifiers is possible, and we plan to include them and show the comparison between our methods and others in revised version.
6. We will compare our method to at least one Shapley-value-based approach, acknowledging that Shapley methods are standard baselines in feature attribution.