Thank you for your detailed review and constructive feedback. We also appreciate your recognition of our method's key strengths: the novel use of normalizing flows to model counterfactual desiderata (including plausibility and actionability) and our demonstration that the proposed method achieves significant computational efficiency without compromising counterfactual quality.

We address your concerns below:

Presentation and Method Clarity

We thank the reviewer for raising this important point. In DiCoFlex, counterfactuals are generated by sampling multiple noise vectors from the base distribution (assumed to be Gaussian) that are further processed using the inverse transformation of the normalizing flow to transform them to counterfactuals. Because we utilize the conditioned variant of the normalizing flow, this process is guided by the initial example $x$, the target class $y'$, and the parameters $p$ and $m$ that control sparsity and controlability. By modeling $p(x'|x,y,p,m')$ as a distribution rather than a point estimate, sampling multiple times naturally produces diverse counterfactuals. The normalizing flow learns the manifold of counterfactuals, and different samples explore different regions of this manifold. That diversity obtained in such a process was empirically confirmed in the evaluation using the Hypervolume metric in our experiments. The role of the score function given by Eq. (6) was just to show how the constraints representing the validity, proximity, and plausibility can be aggregated into one criterion (see more detailed explanation provided to reviewer TN7b). We will revise Section 3 to make the generation process and the role of the score function more explicit.

Furthermore, we provide theoretical analysis of diversity guarantees for our method (See the answer for reviewer TN7b).

Training Details and Hyperparameters

We conducted extensive hyperparameter exploration across the following ranges:

• K: $8, 16, 32$

• alpha: $1, 10, 1000$

• MAF hidden features: $16, 32, 64$

• MAF hidden layers: $2, 5$

• Training configuration: 1000 epochs with early stopping patience of 300 epochs

Training times are reported in a table below:

The final hyperparameters were selected based on validation performance across multiple datasets, balancing counterfactual quality metrics (validity, proximity, plausibility) with computational efficiency. We will include a detailed hyperparameter sensitivity analysis in the final version, showing how these choices affect performance across different datasets and metrics.

Technical Concerns

$L_p$ Norm Stability: Your concern about gradient instability with $p \in (0,1)$ is understandable but doesn't apply here. No gradients pass through this operation. The $L_p$ norm is used only during the k-NN selection phase to define the ground-truth counterfactuals for training. The sparsity control mechanism operates at the data preparation level, not during gradient-based optimization of the flow model.

Related Work: Thank you for highlighting the relevant papers [1,2]. Our key distinction is modeling the conditional probability distribution $p(x'|x,y')$ where both original instance x and target class y' condition the generation. In contrast, the cited works use flows to model $p(x'|y')$ and assess plausibility. Furthermore those methods have limitations regarding processing categorical variables, therefore it is impossible to include them to our benchmark. We will add detailed comparisons to our related work section.

Experimental Questions

Smaller Datasets

We evaluated our method on the German Credit dataset (1,000 samples) as requested. The results demonstrate that our approach maintains excellent performance on smaller datasets:

Our method achieves perfect validity and maintains high classification probability while delivering the best hypervolume performance and competitive LOF scores. These results confirm that our approach scales effectively to smaller benchmark datasets, addressing the reviewer's concern about generalizability across dataset sizes.

Sparsity Explanation

Thank you for seeking clarification on these statements. We provide the following explanation. The statements in Section 3.4 refer to the natural behavior of our k-NN selection mechanism:

Why categorical sparsity is implicit: The key insight lies in the distance metric behavior for different variable types. For categorical variables encoded as one-hot vectors, consider two binary vectors $x=(x_1,x_2)$ and $y=(y_1,y_2)$. The p-distance $d_p(x,y)$ equals 0 when $x_1=y_1$ and equals 2 otherwise, regardless of the chosen p value. This discrete jump property means that to transition between different categories (changing only one categorical variable), the optimization must traverse a fixed distance of 2, naturally limiting categorical changes.

Why continuous features need explicit control: In contrast, continuous variables allow for gradual changes with smaller distances. The choice of p-norm becomes crucial here, as smaller p-values promote sparse changes by concentrating modifications on fewer coordinates. We acknowledge that alternative approaches, like dedicated preprocessing normalization, exist but are beyond the scope of this work.

Limitations and Societal Impact

The reviewer raised an important point about model extraction risks [3,4]. We will expand our limitations section to address potential negative societal impacts, including risk of model extraction through diverse counterfactual explanations.

We believe that these revisions will significantly improve the paper's clarity and completeness while addressing your technical concerns. Thank you for helping us to strengthen this work.

Thank you for your thoughtful and constructive review. We appreciate your recognition of our method's strengths, particularly the model-agnostic approach, dynamic constraint control, and computational efficiency. We address your concerns below:

DiCoFlex vs. s-DiCoFlex Clarification

You are correct that DiCoFlex and s-DiCoFlex represent the same model. s-DiCoFlex is simply DiCoFlex evaluated with p=0.01 during inference to demonstrate sparsity control. Users can dynamically adjust p at inference time. We acknowledge the separate naming is misleading and will revise the paper to clarify that these are configurations of a single model rather than distinct methods.

Sensitivity Analysis Over p Values

We agree that a comprehensive sensitivity analysis strengthens our claims about flexible constraint handling. The table below presents results for different p values ($p \in \{{0.01, 0.08, 0.25, 1.0, 2.0}\}$) for the Adult dataset, demonstrating the smooth control our method provides over the sparsity-proximity trade-off.

The results reveal several key insights: Validity remains consistently high across all p values, confirming robust constraint satisfaction regardless of perturbation level. Clear proximity-sparsity trade-off: Lower p values achieve superior proximity and plausibility, while higher p values yield better categorical sparsity. Smooth transitions: All metrics change gradually across p values, enabling practitioners to fine-tune the balance between competing objectives based on their specific requirements.

This sensitivity analysis validates our method's ability to provide flexible constraint handling, allowing users to systematically navigate the multi-objective trade-off space by simply adjusting the perturbation parameter p.

We will include this analysis to our work in the final version.

Statistical Uncertainty Reporting

We agree that statistical uncertainty measures enhance evaluation rigor. We calculate standard deviations across sets of diverse counterfactual explanations for each instance. We report standard deviation values for all individual-level metrics below. We will add standard deviation values to our main experiment results in Table 1 in the final version.

Note that we exclude Hypervolume from uncertainty reporting as this is a set-level metric that computes a single value per entire counterfactual set, making standard deviation calculations inappropriate.

Sparsity Explanation Correction

Thank you for catching this error. The statement "p = 0.01 encourages lower sparsity" should read "p = 0.01 encourages higher sparsity" (fewer modified features). Lower p values in the Lp norm indeed encourage sparser solutions by penalizing modifications more heavily. We will correct this throughout the paper.

Context Vector Construction Details

The conditioning vector concatenates: $x (original instance), y (target class), p (sparsity parameter), m (actionability mask)$ . This vector is passed to each MAF layer through affine coupling transformations. This way we condition the estimated probability on the conditioning vector as it is described in [1]. We will add a detailed description of this architecture in the Appendix.

[1] Papamakarios, George, Theo Pavlakou, and Iain Murray. "Masked autoregressive flow for density estimation." Advances in neural information processing systems 30 (2017).

Limitations

Since DiCoFlex learns the counterfactual distribution for a specific classifier, classifier changes may invalidate some generated counterfactuals. However, generating multiple diverse counterfactuals allows filtering to retain valid ones, providing natural resilience without retraining.

We believe these revisions will significantly strengthen the paper while maintaining its core contributions. Thank you again for your valuable feedback.

Thank you for your thoughtful review and constructive feedback. We appreciate your recognition of DiCoFlex's clear motivation, thorough experimental evaluation, and computational efficiency advantages.

We address your main concerns below:

Clarification of Score Function Design (Eq. 6)

We agree that introducing the score function given by equation (6) may be confusing and will clarify it in the revised version. This function is not used for training and evaluation. Our intention was to show how the constraints representing the validity, proximity, and plausibility can be aggregated into one criterion and to show that direct optimization of (6) can provide only a single counterfactual explanation, so it does not solve the main goal of this paper, which was focused on generating multiple counterfactuals efficiently.

The connection between validity, proximity, and plausibility and their calculation in eq. (6) is as follows: logprob $\log p(y'|x')$ component represents the correct classification of conterfactual $x'$ to the desired class $y'$ (validity constraint). $d(x',x)$ represents the distance between the original example $x$, and the counterfactual $x'$ (with minus sign in equation (6), represents proximity constraint). $p_{data}(x')$ represents the density function for the data distribution (likelihood function representing plausibility constraint).

As mentioned earlier, we use the score function only for motivation; direct optimization requires access to data distribution (can be modelled with flow) and leads to a counterfactual point, not the distribution of counterfactuals. Therefore, we propose to construct the conditional distribution for generating counterfactual $q(x'|x,y',d)$ given by eq. (7). The generated counterfactuals with this distribution satisfy constraints: Validity, because we select examples from a valid class $y'$. Proximity, because we take the nearest neighbours in class $y'$. Plausibility, because we take the real examples that exist in the data distribution.

Besides, they also satisfy: Sparsity, because various norms can be used for distance calculation. Actionability, because selected attributes can be masked and remain unchanged.

The unparametrized KNN-based model $q(x'|x,y',d)$ is further used to train a parametrized flow-based model $p_{\theta}(x'|x, y)$ by minimizing KLD given by eq. (1) that is equivalent to minimizing $Q$ given eq. (4). The technical implementation is described by the Algorithm in Appendix D (Supplementary Material).

Integration of Score Function and Constraints in Training

Eq. (4) is the criterion used during training and is based on KLD between KNN-based model $q(x'|x,y',d)$ and flow-based model $p_{\theta}(x'|x, y',p,m)$ that is used to train parameters of the flow-based model. The technical procedure for training is described by the Algorithm in Appendix D. In each iteration, we select an example $x$ from the training data and sample a class $y'$. We sample the sparsity level $p$ and mask the candidate $m$. Next, we sample $x'$, one of the $K$ nearest neighbours of $x$ in class $y'$. Note that neighbors are calculated using the distance given by the eq. (9) based on the current $p$ and $m$ values. Finally, the parameters of the flow-based model $p_{\theta}(x'|x, y',p,m)$ are updated with a gradient-based procedure by minimizing $-\log(p_{\theta}(x'|x, y',p,m))$ (based on eq. (4)). We highlight, that we use condition variant of flow-based model, that takes $x$, $y'$ and sparsity parameter $s$ and actionability parametr represetned by mask $m$. So the flow-based model is trained by optimizing Eq. (4) using the neighbours calculated with the distance given by Eq. (9) with $m$ and $p$ sampled in each iteration. Because we sample $m$ and $p$ in each iteration during the training flow-based model, while generating counterfactuals, we can use this model to control sparsity and actionability via simple conditioning on $m$ and $p$.

Theoretical Considerations

We disagree that our approach is largely intuition-driven and lacks rigorous theoretical guarantees. We aim to construct a distribution $q(x'|x,y',d)$ for which the samples represent counterfactuals that satisfy validity (examples from correct class y'), proximity (the closest examples, because we use KNN), plausibility (we select examples that are in data distribution), sparsity (we control the norm with parameter $p$), and actionability (we selected only attributes that can be changed using masking procedure). In lines 175-181, we justify not taking $q(x'|x,y',d)$ directly as a sampling model. Finally, we formulate the theoretical problem of finding a parametrized probabilistic model as a surrogate model in eq. (1), and show how to solve it efficiently.

Additionally, we provide a theorem for theoretical guarantees of diversity preservation for our method that will be included into the final version of our work:

Lemma: Upper bound on difference of expected values between two distributions

Let $p, q$ be two distributions over the same finite set $X$ and the total variation distance $TV(p,q) = \frac{1}{2}\sum_x |p(x) - q(x)|$. If for some function $f: X \to \mathbb{R}$, then:

$\left|\sum_x f(x)(p(x) - q(x))\right| \leq \max_x |f(x)| \cdot \sum_x |p(x) - q(x)| = 2 \cdot \max_x |f(x)| \cdot TV(p,q)$

Theorem: Diversity Preservation in DiCoFlex

Statement

Let $p_\omega(x'|x,y',p,m)$ be the learned conditional normalizing flow and $\hat{q}(x'|x,y',d_{p,m})$ be the empirical distribution defined in Eq. (7). Under assumptions (A1)-(A4), for n samples $x_i', i = 1, 2, \ldots, n$ drawn from $p_\omega$, the expected diversity satisfies:

$$\mathbb{E}\_{p\_\omega}[D(\{x'\_i\}\_{i=1}^n)] \geq \sigma\_{\min}(1 - K^{1-n}) - \Delta\sqrt{2\epsilon}$$

Assumptions

• (A1) The neighborhood $N(x, y', d_{p,m}, K)$ contains $K$ distinct points with minimum pairwise distance $\sigma_{\min} > 0$. Specifically, for any $x_i, x_j \in N(x, y', d_{p,m}, K)$ with $i \neq j$,

$$|x_i - x_j| \geq \sigma_{\min}$$

• (A2) The normalizing flow is trained sufficiently close to the empirical distribution, satisfying:

$$D_{KL}(\hat{q}(x'|x,y',d_{p,m}) \| p_\omega(x'|x,y',p,m)) \leq \epsilon$$

for all $(x,y',p,m)$ in the training distribution.

• (A3) The support of $p_\omega(\cdot\mid x,y',p,m)$ lies in a set of finite diameter $\Delta$, i.e.

$$|x' - x''| \leq \Delta \quad \Longrightarrow \quad 0 \leq D(\{x'\_i\}) = \min\_{i \neq j} |x'_i - x'_j| \leq \Delta.$$

Proof

Step 1: Empirical Diversity Consider the empirical distribution $\hat{q}$. When sampling $n$ points uniformly at random (with replacement) from $K$ distinct points, the expected diversity is at least the product of minimum pairwise distance $\sigma_{\min}$ and the probability of selecting at least two distinct points $P_n(K)$:

$$\mathbb{E}\_{\hat{q}}[D(\{x'\_i\}_{i=1}^n)] \geq \sigma\_{\min} \cdot P_n(K) = \sigma\_{\min}(1 - K^{1-n})$$

Step 2: Total Variation Distance via Pinsker’s Inequality Pinsker's inequality relates the KL-divergence and total variation distance:

$$||p\_\omega(\cdot|x,y',p,m) - \hat{q}(\cdot|x,y',d\_{p,m})||\_{TV} \leq \sqrt{\frac{\epsilon}{2}}$$

Step 3: Bounding the Difference in Diversity The diversity measure, defined by pairwise distances, is bounded above by $\Delta$. Thus, using Lemma 1 we obtain:

$$|\mathbb{E}\_{p\_\omega}[D] - \mathbb{E}\_{\hat{q}}[D]| \leq 2 \cdot \Delta \cdot ||p\_\omega - \hat{q}||\_{TV} \leq 2 \Delta\sqrt{\frac{\epsilon}{2}} = \Delta\sqrt{2\epsilon}$$

Step 4: Combining the Results From Step 3, we have an absolute value bound. Since we specifically seek a lower bound for $\mathbb{E}\_{p_\omega}[D]$, we remove the absolute value by explicitly considering the worst-case scenario for diversity under $p\_\omega$ relative to $\hat{q}$:

$$\mathbb{E}\_{p_\omega}[D(\{x'_i\}\_{i=1}^n)] \geq \mathbb{E}\_{\hat{q}}[D(\{x'_i\}\_{i=1}^n)] - \Delta\sqrt{2\epsilon}$$

Substituting Step 1 into the above, we obtain the desired explicit lower bound:

$$\mathbb{E}\_{p\_\omega}[D(\{x'_i\}\_{i=1}^n)] \geq \sigma\_{\min}(1 - K^{1-n}) - \Delta\sqrt{2\epsilon}$$

Remarks

1. Diversity Guarantee: For large $K$, $1 - K^{1-n}\approx 1$, ensuring strong diversity. Thus the bound

$$\mathbb{E}\_{p\_\omega}[D(\{x'_i\}\_{i=1}^n)] \geq \sigma\_{\min}(1 - K^{1-n}) - \Delta\sqrt{2\epsilon}$$

guarantees that $p\_\omega$ recovers nearly the full empirical diversity except for the training error.

2. Impact of Training: As training improves ($\epsilon \rightarrow 0$), the perturbation term $\Delta\sqrt{2\epsilon}$ vanishes, and the learned distribution’s diversity converges to the empirical distribution’s diversity.