We appreciate the opportunity to clarify your concerns about our work.

W1&Q1 Relationship to PU-Learning

1. Setup: Shared Foundations with PU-Learning

• We acknowledge that our problem shares a foundational setup with PU-learning under the SAR assumption [1,2]. The core challenge of distinguishing "true negatives" from "unobserved positives" is equivalent.
• However, the work in [3] (StableDR) focuses on explicit feedback (e.g., ratings), where "missing data" is less ambiguous. In contrast, our setting involves implicit feedback, where unobserved interactions (no click) are inherently ambiguous—they could stem from non-exposure or disinterest. This ambiguity makes our problem more complex.

2. Method: Divergences from Traditional PU-Learning

• First, our method diverges from PU-learning from the very beginning of problem formulation.
  • Traditional PU-learning [1,2] treats unlabeled samples as a single "unlabeled set" and focuses on weighting or thresholding to identify reliable negatives.
  • We introduce a causal formulation with "exposure" as a treatment variable, stratifying unlabeled samples into three groups (HE, HU, UN) based on estimated exposure confidence. This stratification explicitly models why an interaction is unobserved (high exposure → likely disinterest; low exposure → likely unobserved), which is absent in [1,2].
• Second, our method has a completely novel loss function design.
  • Traditional PU-learning methods such as [1,2] primarily use pointwise losses (e.g., cross-entropy) with reweighting to handle class imbalance.
  • We combine pointwise losses for absolute relavance prediction and pairwise losses to model relative relevance, critical for recommendation ranking tasks. Additionally, we introduce IPM regularization which is to align distributions across exposure states, addressing MNAR-induced shifts ignored in [1,2].

3. Uniqueness of Implicit Feedback in Recommendation Systems

• Recommendation systems aim to rank items by relevance, not just classify them. Thus, pairwise losses (modeling relative relevance between HE, UN, HU) are more aligned with task objectives than the classification-focused losses in [1,2].

W2&Q1 Comparing Generalization Bound of PU learning

Comparison with generalization bound in [4]

• First, we have carefully read the two PU learning-related literatures cited by the reviewers and found that [4] focuses on novel class detection under distribution shift, with its generalization bound constraining the error rate of novel class detection on the target distribution.
• Meanwhile, our work is dedicated to implicit feedback prediction in recommendation systems, and our generalization bound constrains the deviation between the estimated relevance and the true relevance on the whole population.
• The target distribution of [4] is equivalent to the distribution of unlabeled group in our setting. However, the error targets we consider includes both labeled or unlabeled groups, making the two generalization bounds quantitatively incomparable.

Comparison with Generalization Bound in [1]

• The generalization bound in [1] is as follows: $R(\hat y_{\hat R}\mid y)\leq \hat R(\hat y_{\hat R}\mid e,s)+\sqrt{\frac{\delta^2_{\max}\ln \frac{|\mathcal{H}|}{\eta}}{2n}}.$ As a classical generalization bound for empirical risk minimization (ERM), it only depends on the size of the hypothesis space $|\mathcal{H}|$ and the total sample size $n$, which measn that it includes all the error of MNAR and positive unlabelled machenism in the single term.
• Meanwhile, our generalization bound explicitly model four types of error sources:

1. Empirical risk for DP and HE groups,
2. Distribution shift between exposed ($O=1$) and unexposed ($O=0$) groups with high confidence ($R=1$),
3. Distribution shift between high confidence ($R=1$) and low confidence ($R=0$) groups,
4. Irreducible variance of potential outcome in the exposed ($O=1$) and unexposed ($O=0$) groups.

• Our generalization bounds provide a more detailed analysis of the sources of error, while [1] lumps all sources of error together and bounds them using the Hoeffding inequality.
• Intuitively, in the following situations, our bound will have a more significant advantage compared to [1]:

1. When the MNAR problem is severe, the IPM term in our bound directly quantifies the distribution difference between the exposed and unexposed groups, making it more suitable for the popularity bias scenarios commonly seen in recommendation systems.
2. When the proportion of unexposed samples is high, our bound can sense this proportion by weighting the distribution shift of low-confidence samples (UN group), which is more in line with the cold start scenario in the recommendation system.
3. When the variance of potential outcomes is small, the counterfactual prediction has high certainty, and the upper bound of our bound will be significantly reduced, which is more in line with the scenario in the recommendation system where users have similar strong preferences for specific items. In addition, empirically we also compared our method with some PU learning methods under the SAR assumption on the Coat and Yahoo datasets. The following experimental results demonstrate the effectiveness of our method.

Bekker [1]

PGPU

PUe

W3&Q2 Motivation of Algorithmic Choices

• Inclusion of pairwise losses. We are motivated by the Bayesian Personalized Ranking (BPR) [5] from implicit feedback community. When the label for each sample is not exactly accurate, while we are more confident about the relationship between two samples, a pair-wise loss is suitable. The effectiveness of inclusion of pairwise losses can also be validated by the following ablation study.

• Definition of D_+. We actually have included the definition in line 219 of the original manuscript. We use D_+ to represent the positive group in each of our three pairwise losses.

• Sigmoid in Pair-wise Losses. We follow the convention of BPR loss [5] to use a sigmoid function to normalize the difference between pairs. Intuitively, if the sigmoid was omitted, some pairs with significant difference would dominate the pairwise losses.

• Hyperparameters Tuning. It is easy to find that for a given dataset, there is a 1-1 correspondence between $\alpha$ and $D$, and also for $\beta$ and $D'$. Besides these two hyperparameters, there is also a kernel parameter $q$ and the loss weights $\lambda_{point},\lambda_{pair}$. The first three hyperparameters are tuned in the first stage (stratification), while the last two are tuned in the second stage (training), which will significantly reduce the searching space.

W4 RBF Kernel and Other Type of Datasets

We acknowledge that the RBF kernel in the paper is a straightforward and naive approach. However, in fact, our method can incorporate more efficient kernel estimation methods easily. We added the experiments using random Fourier feature and Nystrom approximation, which shows that the performance is retained while the runtime reduces to one-fourth of the original SVDD.

In addition, based on this, we transfer to an image-type dataset and compared the breast cancer detection performance with some PU learning methods. The experimental results are as follows.

[5] Rendle, S., et al. 2009. BPR: Bayesian personalized ranking from implicit feedback. In Proceedings of the Twenty-Fifth Conference on Uncertainty in Artificial Intelligence (UAI '09).

[6] He, F., et al. "Instance-Dependent PU Learning by Bayesian Optimal Relabeling. arXiv, Aug 2018." URL https://arxiv. org/abs (1808).

[7] Wang, Xutao, et al. "Pue: Biased positive-unlabeled learning enhancement by causal inference." Advances in Neural Information Processing Systems 36 (2023): 19783-19798.

Again, we sincerely thank you for the insightful comments.

During the rebuttal period, we have validated the effectiveness of our method compared with some state-of-the-art PU learning baselines under SAR assumption, and we empirically showed that our method performs even better with other kernel estimation methods such as random Fourier feature and Nystrom approximation. Furthermore, we added another experiment on a public image-type dataset for breast-cancer detection. We have conducted the experiment on the Stanford Sentiment Treebank (SST-2) where the goal is to classify positive and negative reviews. Due to the large scale of the dataset and the backbone, we have to present these results now.

Specifically, the original SST-2 dataset has 68% positive samples and 32% negative samples. We conduct experiments under both SCAR and SAR assumptions. To simulate the SCAR setting, we randomly flip some positive labels to negative according to a fixed probability $p$. We choose various values for $p$, i.e., 0.1, 0.3, 0.5, and the experimental results are as follows.

Our method stably outperforms the PU learning baselines under various levels of unlabeled probability.

To simulate the PU setting under SAR assumption, we use a pre-trained RoBERTa to get the representations with $d$ dimensions for all sample inputs $x_i\in\mathbb{R}^d$. And for the training and development set, we flip some positive labels to negative according to the following probability:

$p^{flip}_i=\sigma(W\cdot x_i),$ where $W\sim \mathcal{N}(0, 1)^d$ is a random vector applied for all samples. Note that the labels of test set is not flipped since we want a concrete evaluation of classification performance. The experimental results are as follows.

Image-type and text-type datasets provide a more challenging setup compared to interaction matrix. However, the proposed method stably outperforms the previous PU learning baselines, showing the transferability of our method to various type of datasets.

We hope these extra experimental results help to solve your concerns on transferability. And we would be happy to provide further clarification or address any additional questions you may have regarding our manuscript. Please do not hesitate to let us know if there is anything else you wish to discuss.

Hence the main novelties I see are in the specific way the proposed method identifies and uses these reliable examples, in formal results (generalization bound) and empirical performance. These may be reasonable contributions, but the paper is currently written in a way that obfuscates contributions and emphasizes aspects that are not technical novelties but rather are choices of notation.

In addition to the above two methodological novelties, we consider our derived formal results (in particular the generalization bound) as our main technical contributions. We also agree with the reviewer that it's necessary to further compare our generalization bound to those in the PU-learning literature under SAR (see details below).

I appreciate the comparison of generalization bounds included in the rebuttal and think a summary of these comparisons should appear in the paper.

Thank you so much - we will definitely add a summary of these comparisons included in the rebuttal to our final version. This helps to clarify the theoretical importance of our findings.

Yet I disagree with some of the claims made in the rebuttal. It is claimed that the generalization bound in [4] cannot be compared with the one in this work, since it bounds the error over the target/unlabeled distribution and not the entire population. But this is a minor technical point, it is easy to add an error bound term for the labeled population with standard empirical risk minimization results.

We highly appreciate your comments - as suggested by the reviewer, we admit it's possible to make fair comparison by adding an error bound term for the labeled population with standard empirical risk minimization results. We would like to share with you the following interesting findings. Making such comparison helps clarify the key differences with our work. Here, we elaborate on the key distinctions and connections.

It is easy to add an error bound term for the labeled population with standard empirical risk minimization results. The main difference seems to be in the divergences used by the two results, and some mention of these differences would be useful.

We can frame the comparison using the law of total probability for the overall risk, where the total population is a mix of the source ($Y=0$) and target ($Y=1$) distributions:

$Risk_{total}=P(Y=1)\cdot Risk_{Y=1}+P(Y=0)\cdot Risk_{Y=0}.$ Error Bound on the Source Population

For the source population ($P_S$), where we have precise labels (i.e., known non-novel instances), the error of a classifier $\hat h$ is the false positive rate. This can be bounded using standard Empirical Risk Minimization (ERM) theory. The risk is bounded by the empirical risk on the sample $P_S$ plus terms for model complexity and statistical confidence:

where $\hat{E}_{S_S}$ is the empirical error on the source set $S_S$ and $\Re$ is the Rademacher complexity.

Error Bound on the Target Population

For the target population ($P_T$), which contains a mix of the shifted nominal distribution ($P_{T,0}$) and the novel category ($P_{T,1}$), Theorem 4.3 from [4] provides a direct bound on the classification risk. $\text { Error Bound on } Y=1: \quad R_T^{l_{01}}(\hat h) \leq R_T^{l_{01}}(h^\*)+4 \epsilon_{\text {shift }}+2(\beta-\beta(h^\*))+\Re_{n_S,P_S}(\mathcal H)+\Re_{n_T,P_T}(\mathcal H)+\sqrt{2\ln(1/\delta)}[n_S^{-\frac{1}{2}}+n_T^{-\frac{1}{2}}].$

The key term here, $\epsilon_{shift}$, is a bound on the constrained H-divergence, $d_{\mathcal{H},\beta}(P_S||P_{T,0})$, which measures the maximum discrepancy on "rare" events.

Overall Bound for the Framework in [4]

Combining these, the overall risk for the method in [4] across the entire mixed population (with mixture proportions $\pi_S$ and $\pi_T$) can be bounded as:

$Risk_{total}^{[4]} \leq \pi_S ( \hat E_{S_S}[\hat h(x)] + C_1) + \pi_T ( R_T^{l_{01}}(h^\*) +2(\beta-\beta (h^\*))+C_2),$ where $C_1=2 \Re_{n_S, P_S}(\mathcal H)+\sqrt{\frac{\ln (1 / \delta)}{2 n_S}}$, and $C_2=\Re_{n_S,P_S}(\mathcal H)+\Re_{n_T,P_T}(\mathcal H)+\sqrt{2\ln(1/\delta)}[n_S^{-\frac{1}{2}}+n_T^{-\frac{1}{2}}].$

Our Overall Bound and Comparison

Our paper directly provides a single, unified bound on the mean squared error for the potential outcome $Y(1)$ over the entire population, which is as follows:

$Risk_{total}^{ours} \leq \underbrace{2E_{x\mid r, o}[(h(\Phi(x))-Y(1))^2\mid R=1, O=1]}_{\text {factual loss of the DP and HE groups}}$

$+ 2B_\Phi \underbrace{IPM_G(p_\Phi^{R=1,O=0}, p_\Phi^{R=1,O=1})}_{\text {distribution shift on O given R=1}}$

$+4 P( R=0 ) B_\Phi \underbrace{IPM_G (p_\Phi^{R=0}, p_\Phi^{R=1})}_{\text { distribution shift on R}}$

$-8 \underbrace{(\min (\{\mathbb V(Y(1)|O=0),\mathbb V(Y(1)|O=1)\}))^2}_{\text { minimal variance of potential outcomes}}.$

Thank you for your constructive remarks and valuable feedback! Below, we address your questions and indicate the changes we’ve made thanks to your suggestion (I am so sorry for my late reply, due to serious medical/family emergency).

My interpretation of the paper is that the use of counterfactual language is just a choice of notation, and it is independent of any technical contribution. I do not object this choice of notation (it makes the no-unobesrved-confounding assumption explicit), but once the assumptions are made the problem is formally equivalent to a standard PU-learning problem.

For problem setup and notations, we agree that the formulated counterfactual problem is formally equivalent to a standard PU-learning problem under MAR. In fact, the studied implicit feedback problem in recommender system community is essentially a standard PU-learning problem under MAR in machine learning community.

As confirmed in the rebuttal, this notation is also not novel, hence I do not understand why it is claimed as a contribution, moreover a major one. Other reviewers pointed out the use of counterfactual language as a signifiant strength and novelty of the paper, so I'd be interested to hear their updated view after reading the discussion here.

We have reservations about the novelty of formulating the implicit feedback as a counterfactual inference under missing treatment problem. Please kindly note that we agree the core contribution is not notation, but technically motivates a new counterfactual inference method for addressing the implicit feedback problem with strong theoretical guarantee (see parts 2 and 3 for our detailed discussion). This should be seen as a technical contribution to the causal community, observing that pervious counterfactual imputation error bound only limited to the standard counterfactual inference with fully observed treatment status.

To my best understanding, works on identifying reliable negatives, which I referenced in the review, also have the same basic goal and intuition as this work. The authors did not comment on this point in their rebuttal.

We fully agree that working on identifying reliable negatives (in our first stage) in standard in PU-learning. The main technical (including both methodological and theoretically) contributions instead are in our second stage (see below).

• For facilitate checking, we recap our stratification table in below with proper comments. The $\textcolor{red}{red}$ font indicates that the value is imputed.

After identifying reliable negatives in the first stage, our novelties are in the second stage that:

• Methodologically, we propose a representation learning-based counterfactual inference method to address the missing treatment issue, motivated by BNN, BLR [1], we use $\mathcal L_{IPM}= IPM_G(p_\Phi^{R=0}, p_\Phi^{R=1}) + IPM_G(p_\Phi^{R=1, O=0}, p_\Phi^{R=1, O=1})$ to:
  • balance the distribution shifts between unreliable treatments ($R=0$, including UN) and reliable treatments ($R=1$, including DP, HE, and HU);
  • among reliable treatments ($R=1$, including DP, HE, and HU), balance the distribution shifts between treated ($O=1$, including DP and HE) and untreated units ($O=0$, including HU).

[1] Uri Shalit, et al. Estimating individual treatment effect: generalization bounds and algorithms, ICML 2017.

• Methodologically, we propose a contrastive learning loss to estimate the counterfactuals under missing treatment, by using $\mathcal L_{pair} (h(\Phi(X_{DP})), h(\Phi(X_{HE})))$+$\mathcal L_{pair} (h(\Phi(X_{HU})), h(\Phi(X_{UN})))$+$\mathcal L_{pair} (h(\Phi(X_{UN})), h(\Phi(X_{HE})))$, for reasons below:
  • As shown in the above Table, $Y(1)=1$ in DP while $Y(1)=0$ in HE, so we set (DP, HE) as positive-negative pairs using $\mathcal L_{pair} (h(\Phi(X_{DP})), h(\Phi(X_{HE})))$;
  • Importantly, $P(O=1)$ in UN is higher than in HU, i.e., UN group has higher exposed rate than HU, so users give negative response in UN maybe due to they are less interest to the items compared with HU, i.e., $Y(1)$ in UN has less chance to be 1 compared with HU, leading to $\mathcal L_{pair} (h(\Phi(X_{HU})), h(\Phi(X_{UN})))$;
  • Following a similar argument, we have $\mathcal L_{pair} (h(\Phi(X_{UN})), h(\Phi(X_{HE})))$.

Please kindly note that both (i) representation learning and (ii) contrastive learning are in our second stage, that is, after identifying reliable negatives.

The key distinction is the divergence used to control for distribution shift.

Our bound is more insightful and applicable for the specific problem of implicit feedback for recommender systems. Our bound provides a superior analytical framework in the following key scenarios:

1. When MNAR is the Dominant Problem.

The core challenge in implicit feedback is often selection bias from MNAR scenario. This is a global distribution shift. The IPM in our bound is designed to capture such global geometric shifts, whereas the constrained H-divergence in [4] focuses on rare events and is less sensitive to global shifts. Our bound's specific term for MNAR bias (IPM distance) directly account for this primary challenge.

2. When Granular Diagnosis of Bias is Required.

Our bound's decomposition allows for fine-grained risk analysis, especially for the distribution shifts. The single $\epsilon_{shift}$ term in [4]'s bound identifies that a distribution shift exists but does not offer insight into its nature or source, making it less of a diagnostic tool.

3. When unlabeled samples have strong heterogeneity.

In implicit feedback, non-clicks ($Y = 0$) may be due to "non-exposure" (HU group) or "irrelevance" (HE group), and the meanings of these two for prediction are completely opposite. Our bound constrains the errors separately through stratification (HE/HU/UN) to avoid misclassifying "potentially positive samples that are not exposed" as negative samples; the bound in [4] treats all unlabeled data as homogeneous, cannot distinguish between "non-exposure" and "irrelevance", and will overestimate the negative sample error.

In the ideal scenario when the divergences used by the two results are all zeros, i.e., $\epsilon_{shift}=0$, $IPM_G(p_\Phi^{R=1,O=0}, p_\Phi^{R=1,O=1})=IPM_G(p_\Phi^{R=0}, p_\Phi^{R=1})=0$, the bounds are reduced to the following forms:

Risk_{total}^{[4]} \leq \pi_S ( \hat E_{S_S}[\hat h(x)] + C_1) + \pi_T ( R_T^{l_{01}}(h^\*) +2(\beta-\beta\(h^\*))+C_2),

$Risk_{total}^{ours} \leq \underbrace{2E_{x\mid r, o}[(h(\Phi(x))-Y(1))^2\mid R=1, O=1]}_{\text {factual loss of the DP and HE groups}}$

-8 \underbrace{(\min\{\mathbb{V}(Y(1)|O=0),\mathbb{V}(Y(1)|O=1)})^2}_{\text { minimal variance of potential outcomes}}.

• Bound [4] measures the gap between the learned model $\hat h$ and the ideal model $h^*$. The goal is to see how close you can get to the absolute best performance possible.

• Our Bound measures the gap between the model's performance on the entire population and its performance on the observed subset (the DP and HE groups). Its structure is typical of causal inference and semi-supervised learning. The goal is to see how well you can generalize from what you know to what you don't.

• Bound [4] has an irreducible error term $R_T^{l_{01}}(h^\*)$, the Bayes error rate. This is a hard floor on performance; no model can do better.

• Our Bound has a "bonus" term, $\min( {\mathbb{V}(Y(1)|O=0),\mathbb{V}(Y(1)|O=1)})$, which is related to the variance of potential outcomes. This is not a performance floor. Instead, it indicates that if user preferences are deterministic (i.e., variance is low), the learning problem itself becomes easier, and the bound tightens to reflect this. This concept is native to potential outcome frameworks.

Please let us know if we have properly addressed your questions and we are more than happy to discuss more!

Thank you for these detailed and constructive suggestions. Here are our responses.

W1 Hyperparameter Analysis

We have actually conducted the sensitivity analysis for the DP sample ratio hyperparameter and the negative sample ratio hyperparameter as Figure 3 of the submitted manuscript shows. Unfortunately, we made a typo in Fig 3, which may have misled you. The $D$ in Fig 3 should be the corresponding $\alpha$ and $P_{HU}$ should be $\beta$ defined in Section 3.3.

W2 Formatting Issue

Thanks to point out, we have made a typo here, and the correct sentence should be In contrast, when the samples are not assigned the treatment (i.e., $O_{u,i} = 0$, corresponding to UN and HU groups), we treat these unobserved interactions as counterfactual data and employ pairwise loss to model the relative rankings.

W3 Figure 2 Explanation

Thanks for your suggestion. Figure 2 illustrates our causal representation learning framework: input features are mapped via $\Phi$ to a latent space, where pointwise losses (DP/HE) and pairwise losses (UN vs. HE, HU vs. UN) guide training. IPM regularization aligns distributions of $R=0$ and $R=1$ samples, as well as the $R=1, O=0$ and $R=1, O=1$ groups. We will modify the manuscript so that each part corresponding to the figure is refered.

W4 Counterfactual Baselines

Thanks for your suggestion, we added comparisons with IPS and DR on Coat and Yahoo datasets. Our method outperforms them significantly. We will include these baselines in the final version of this paper.

W5 Generalization Bound Explanation

• The first term is the factual loss based on the true value of $Y(1)$ of the DP and HE groups. It measures the error of potential outcome estimation for those in DP and HE groups.

• The second and third terms are the IPM distance measuring the distribution shift on $O = 1$ and $O = 0$ given $R = 1$ group and distribution shift on $R = 1$ and $R = 0$ groups weighted by the proportion of UN group $P(R = 0)$.

• The last term measures the minimal variance of potential outcomes on $O=0$ and $O=1$ groups, which is independent of the model selection.

• This bound actually guides the loss design for our method. Intuitively, minimizing the L_point and L_pair can control the factual loss, while minimizing the L_IPM is an effective approach to handle the distribution shifts.

W6 Appendix Clarification

We appreciate the reviewer for pointing out this flaws. We should clarify that the Related Work is in Appendix A, the proof for our theory is in Appendix B, the data preprocessing protocal is in Appendix C, the experimental results on KuaiRec dataset is in Appendix D, while Appendix E provides the compute resources we used to conduct all the experiments.

W7 t-SNE Visualization

As suggested by the reviewer, we conduct experiments on a synthetic dataset and visualize the distributions of the DP, HE, HU, and UN groups with t-SNE plots. As expected, we see the distribution of HE and HU groups aligns well with the exposed and unexposed populations, and the UN groups lie in the boundry of the two populations. Unfortunately, due to the rule of NeurIPS 2025, we cannot share the visualization results by images at this stage. We are delighted to include the figures in the final version of this paper.

Q1 Cold-Start Suitability

Thank you for providing another potential benefits of our method. We conduct experiments for cold-start users and compare with existing baselines. We investigate the users with the lowest 5%, 10%, 15% or 20% of interactions, and analysis the drop against normal setting. Performance on cold-start dataset is slightly lower than normal setting but still outperforms baselines.

Q2 Notations

No they are not. In Section 2, we defined $D_1=\{(u,i)|(u,i)\in D,Y_{u,i}=1\}$ and $D_0={(u,i)\mid(u,i)\in D,Y_{u,i}=0}$, which is only based on the observed outcome to devided the dataset. While in Section 3.4, the positive and negative samples refer to two groups present in the pairwise loss, which is based on our stratified dataset. We have included 3 pairwise losses, so the positive and negative samples may refer to (DP, HE), (UN, HE) and (HU, UN) respectively.

Q3 Imputation Errors

Thanks for providing another in-depth analysis perspective for our method. We compare the proposed imputation method with random imputation and K Nearest Neighbor (KNN)-based imputation on Coat and Yahoo dataset. Here we choose K=3 and impute a unlabeled sample as positive if all 3 neighbors are positive, and negative if all 3 neighbors are negative. The experimental results shows that our kernel-based imputation significantly outperforms the other imputation methods. This indicate that the imputation quality has a significant impact on the subsequent prediction results, and our kernel-based imputation is effective.

We really appreciate your recognition of our paper's contributions and your kind words, and we are happy to hear that your concerns have been addressed - thank you so much!

Gaussian kernel introduces significant computational overhead

Thank you for highlighting these critical points. We acknowledge that the kernel-based hypersphere model introduces $O(m^2)$ complexity with $m$ positive samples, which poses challenges for large-scale data. We have explored scalable alternatives as you suggested:

Random Fourier Features (RFF) and Nyström approximation reduce kernel matrix computations from $O(m^2)$ to $O(md)$, where $d$ is the dimension of features and is significantly smaller than $m$. Preliminary tests on the Coat and Yahoo dataset show that RFF performs even better than original SVDD while cutting runtime by xx%. Nyström retain xx% of the original performance while cutting runtime by xx%. We plan to include this extension in the revised manuscript.

We can now achieve to address the scalability concern for this paper by replacing the kernel-based method. Sincerely thank you again for your valuable suggestions, which make our work much more solid.

Thank you so much for your kind suggestions - we will definitely include the above additional experiments on RFF and Nyström approximation in our final version. Many thanks for holding a positive opinion of our paper!

We appreciate your insights and address each point below.

Unconfoundedness Assumption Although the violation of unconfoundedness assumption is arthogonal to our contribution, and the unconfoundedness assumption is common in implicit feedback recommendation litertures [1,2], we agree with the reviewer that it is meaningful to validate the robustness of our proposed method against its violation. To test robustness, we conduct experiments on a synthetic dataset, introducing varying levels of unmeasured confounding. The data generating process is as follows.

• $U\sim \mathcal{N}(\mu,1)^{n_U\times n_I}$,
• $\tilde r_{u,i}=U_{u,i}+\epsilon_{u,i}$,
• $r_{u,i}=round(\min(0,\max(\tilde r_{u,i},5)))$,
• $O_{u,i}\sim Bernoulli(\sigma(U_{u,i}-\mu))$

Here the potential outcome and the missing mechanism is counfounded by an unmeasured variable $U$. We take $\mu=1,2,3$ for different strength of unmeasured confounding. The results show that our proposed method consistently outperforms the baselines and indicates reasonable resilience against stronger confounding. We will include these results in the revised version of the manuscript.

SVDD Complexity and Scalability

The original SVDD step has $O(N^2)$ complexity, but we mitigate this by implementing the Random Fourier Features (RFF)** and Nyström approximation reduce kernel matrix computations from $O(N^2)$ to $O(Nd)$, where $d$ is the dimension of features and is significantly smaller than $N$. Surprisingly, we find that they even outperforms the original SVDD while reduces the runtime by 75%. We will include these results in the revised version of our manuscript, and hope this can address your concern about scalability.

Compare with RFF and Nystroem

Hyperparameter Tuning

Despite our method has 5 hyperparameters in total, but are used in two separate stages, i.e., the percentiles for stratification $\alpha, \beta$ and kernel parameter $q$ at the imputation and stratification stage, while loss weights $\lambda_{point}$ and $\lambda_{pair}$ in the counterfactual representation model training stage. In this way, the searching space is much smaller.

We are happy to hear that your concerns have been addressed. Thank you for your efforts in reviewing our paper!