ReTaSA: A Nonparametric Functional Estimation Approach for Addressing Continuous Target Shift

4. There are some gaps between the theoretical results of Section 4 and the real-world OOD generalization problem the paper seeks to solve: 4.1. Theorem 1 bounds the $L^{2}$ error of the estimated $\rho$. However, $\rho$ is only a means to re-weighting the risk function to adapt to the test domain (as explained on Page 3), and it's not clear to me whether estimating $\rho$ well in $L^{2}$ distance is necessary or sufficient to adapt to the test domain. I think the paper should provide some more concrete connection between estimation of $\rho$ in $L^{2}$ distance and test-domain performance of the new risk minimizer.

Our Response: Thanks for pointing this out. The key reason for the weight estimation is motivated by the loss relationship between training and testing datasets.

With a well-estimated $\rho$, we can recover the importance weight function $\omega$ accurately, which can facilitate us to learn a predictive model with respect to the test-domain distribution. As $\rho$ is a function, we use the $L^2$-distance to measure the estimation accuracy in the functional space.

4.2. It's unclear (to me) how some of the main assumptions in Section 4 relate to the real-world problem being solved; see Major Questions 3. and 4. below.

Our Response: Thanks for the comment. Please see our responses below under Major Questions 3-4.

Weaknesses:

Major

1. The paper should discuss the theoretical computational complexity and practical scalability of ReTaSA, in terms of the source and target sample sizes $n$ and $m$, data dimension $p$, etc. Relatedly, under "Related Work", the paper claims "empirical evidence from ... our experimental studies confirms that KMM is computationally inefficient in categorical and continuous cases," but I couldn't find evidence of this in the paper.

Our Response: Thanks for pointing this out. Our method computationally faster than the KMM method is because our method only needs to solve a linear equation, while KMM needs to optimize the constraint quadratic programming problem. More specifically, the major computational cost in our method is inverting a $n\times n$ matrix, whose computational complexity is generally $\mathcal{O}(n^3)$. The KMM method proposed by Zhang et al. (2013) needs to optimize a quadratically constrained quadratic program with a $n\times n$ matrix. We believe the computational efficiency can be attributed to the fact that our proposed method only requires a one-step linear system solving while the KMM method (as well as the added L2IWE method) needs to loop through iteration in the optimization process. From the empirical aspect, we detail the running time comparisons in Figure S.4 in Section C.2.2 of the supplementary materials.

2. I found several parts of the paper a bit vague or missing some details that would help the reader: 2.1. Page 2, last sentence, "the target shift assumption implies that... (1)": I think it would be helpful to include a few more details on the steps by which the target shift assumption implies Eq. (1). I was eventually able to figure this out (using Bayes' rule), but it interrupted my reading of the paper and took a few minutes. This could easily be avoided by adding another intermediate equality in Eq. (1) without increasing the length of the paper.

Our Response: Thanks for the suggestion. We have provided more details in the derivation in the revision in Equation (1).

2.2. Page 3, just after Eq. (3), "$T$ and $T^{\ast}$ are adjoint operators because $\langle T\phi,\psi\rangle = \langle\phi, T^{\ast}\psi\rangle$": It's not immediately obvious why this is the case. Since this observation is important for the remainder of the paper, please include a more detailed explanation or proof in the main text or indicate where this could be found (e.g., in an appendix).

Our Response: We have refined the definition of inner products and added a few more lines of equations between equation (3) and (4) to validate this fact.

3. Page 3, just after Eq. (4), "where... $\rho(y)$ is an unknown function to be solved": I found this quite confusing because it sounded like $\rho = \omega$. Only later is it explained that $\rho = \omega - 1$. So perhaps this latter fact can be explained a few sentences earlier.

Our Response: Thanks for the suggestion; we have clarified the usage of notation $\rho$ as you suggested and edited the identifiability subsection to reflect this change.

3. First Paragraph of Section 3: $\eta = p_{t}/p_{s}$ is estimated as the ratio of two density estimates. There is a significant body of work on density ratio estimation showing that estimating the ratio by the ratio of two density estimates is often suboptimal, both in theory and in practice (see, e.g., $K17, SSK10$ ). Given this, perhaps the paper should consider using direct density estimation methods for this step, both to improve practical performance and to relax the assumptions (specifically, Assumption 3).

Our Response: Thanks for providing closely connected references. Utilizing direct density estimation to yield an improved rate of consistency will be an interesting future theoretical work to pursue. In the meantime, we believe that assumption 3 in our setting is natural as the equation (1) would not be valid if violated. We have added a remark to explain connections between the assumption 3 and the equation (1).

Furthermore, in the revised paper, we compared two ways of computing the density ratio: directly computing the ratio of two estimated densities and the Relative unconstrained Least-Squares Importance Fitting (RuLSIF) method (Yamada et al. (2011)), which is proposed by the author of [SSK10]. The numerical study shows that these two methods have similar performance. The details of the numerical comparison can be found in Section C.2.3 of the Supplementary Materials.

Thanks to the authors for responding to my comments. For my comment 4.2, I was hoping for was some more detailed discussion along the following lines.

Ultimately, rather than estimating $\omega$, what we care about is something like the test-domain generalization error $E_s[\omega(y) \ell(\hat f(x), y))] - \frac{1}{n} \sum_{i = 1}^n \hat\omega(y_i) \ell(\hat f(x_i), y_i)$ of the weighted risk minimization $\hat f = \arg\min_f \frac{1}{n} \sum_{i = 1}^n \hat\omega(y_i) \ell(f(x_i), y_i).$ One approach to bound this generalization error is by bounding the difference between the true and empirical risks uniformly over the hypothesis class, i.e., $\sup_f \left| \frac{1}{n} \sum_{i = 1}^n \hat\omega(y_i) \ell(f(x_i), y_i) - \mathbb{E}_s[\omega(y) \ell(f(x), y))] \right|.$

For simplicity, let's assume the source sample size $n$ is huge, the hypothesis class is reasonably bounded, and $y$ is bounded almost surely, so that the true and empirical risk functionals are basically identical; i.e., $\sup_f \left| \frac{1}{n} \sum_{i = 1}^n \hat\omega(y_i) \ell(f(x_i), y_i) - \mathbb{E}_s[\hat\omega(y) \ell(f(x), y)] \right| \approx 0$ (this can be bounded more formally by standard techniques). Then, the generalization error mostly comes from the error in estimating $\omega$, namely

Our Response: Thanks for the clarifications. The singular value decomposition (SVD) of general conditional expectation operators is non-trivial compared with its counterpart for matrices. To provide some intuitive explanations, we focus on the exponential family. More specifically, the Gaussian distribution in our response. Suppose we can factorize a joint distribution $p(x,y)$ into $p(x|y)p(y)$, where

In fact, this joint distribution $p(x,y)$ corresponds to the model $x=y+\epsilon$, where the random variable $y\sim\mathrm{Norm}(0,\sigma_0^2)$ is independent of the random Gaussian noise term $\epsilon\sim\mathrm{Norm}(0,\sigma^2)$. Through SVD of the conditional expectation operator $T$ (associated with $p(x|y)$), the singular values have a closed form and are given by (see Makur and Zheng 2016 for more details)

We can see that the decay rate of $\lambda_i$ depends on the ratio $\sigma^2/\sigma_0^2$. Recall the model behind the joint distribution is $x=y+\epsilon$, where $\sigma^2$ is the variance of $\epsilon$ and can be seen as noise; and $\sigma_0^2$ is the variance of $Y$ and can be seen as information. Thus, we can see that if the noise dominates the information (i.e., $\sigma^2\gg\sigma^2_0$), the decay is fast. On the other hand, if the information dominates the noise (i.e., $\sigma_0^2\gg\sigma^2$), the decay is slow. If $\beta$ is given, and the noise-to-information ratio is high, the $\beta$-regularity space becomes smaller. In other words, it becomes more difficult to identify functions in the high-noise setting. The reason of the resulting smaller space is that we require the Fourier coefficients to decay at a fast rate to ensure the series is finite.

Now, back to the general model $p(x,y)$, qualitative speaking, the singular values $\lambda_i$ of the operator $T$ describe how well we know about $x$ based on the information of $y$. If $y$ does not provide much information about $x$, or in other words, $x$ contains a lot of noise, the singular values decay fast, and we say the conditional expectation operator is severely ill-posed, and vice versa.

The $\beta$-regularization spaces, as we have explained in the previous responses, are to characterize both the singular values of the conditional expectation operator $T$ and the Fourier coefficients of the function $\rho(\cdot)$. The appropriate value for $\beta$ also depends on both $T$ and $\rho(\cdot)$. A larger $\beta$ implies a smaller function space for $\rho(\cdot)$ to live in, given the operator $T$. Lastly, under our identifiability assumption, we would like to point out that the $\beta$-regularity space can also be viewed as the range of the operator $(T^\ast T)^{\beta/2}$.

Reference

Makur, A., & Zheng, L. (2016). Polynomial spectral decomposition of conditional expectation operators. In 2016 54th Annual Allerton Conference on Communication, Control, and Computing (Allerton). IEEE.

Our Response: Thanks for your thought-provoking insights. First, we apologize for the confusion caused by our notation and terminology. Our norm is different, unlike the usual $L^2$ norm of functions, defined as $\Vert f\Vert=\sqrt{\int f(x)^2 dx}$, we define the inner product for the Hilbert space as $\langle f(x,y), g(x, y)\rangle = \int f(x,y) g(x, y) p_s(x,y)dxdy=\mathbb{E}_s(f(x,y)g(x,y))$ and norm as $\Vert f\Vert=\langle f, f\rangle^{1/2}$ (see page 3 of our paper). Thus, the norm we use for a function $f(y)$ is defined as $\Vert f\Vert=\sqrt{\mathbb{E}_s(f(y))^2}$, which is exactly the “weighted $L^2$ norm” you mentioned in your comments.

Dear Reviewer SE3s,

We extend our sincere appreciation for your invaluable feedback and constructive suggestions. As we approach the conclusion of the rebuttal/discussion phase, we would like to inquire if there are any remaining questions or concerns that we can address or clarify.

If you find our revision and response satisfactory, could you please consider raising the score to support our efforts? Your support for our work and insightful suggestions for improving its quality have been immensely valuable to us. We are grateful for the time and consideration you have dedicated to our submission.

Thank you once again.

Authors.

We sincerely thank you for your invaluable feedback and guidance throughout the review process. Your insights have been incredibly helpful in improving our paper, and we have implemented all the suggestions you provided accordingly in the revision: We added the discussion on 1) $\beta$-regularity in supplementary Section A.1 and 2) the relationship between the weighted $L^2$ norm and test-domain generalization error in supplementary Section B.1.

Once again, thank you for your time, commitment, and valuable suggestions. We appreciate your dedication to ensuring the highest standards and are grateful for your support.

Comment 1

I found the major issue is insufficient comparison with other prior works. For example, Nguyen et al., (2016) studied the same problem with a slightly different estimation procedure. It would be nice to discuss and compare the method.

While the estimation procedure is straightforward, it would be nice to clarify what is the technical challenge and the novelty of this method.

Nguyen, T. D., Christoffel, M., & Sugiyama, M. (2016, February). Continuous target shift adaptation in supervised learning. In Asian Conference on Machine Learning (pp. 285-300). PMLR.

Our Response

In our revised paper, we added the L2IWE method proposed by Nguyen et al. (2016) for extensive comparisons. Our proposed ReTaSA method still outperforms the L2IWE method consistently. The KMM and the L2IWE methods are similar in spirit as they both try to estimate the importance weight function through distribution matching; their difference is that they choose different function spaces to search for the importance weight function. The novelty of our method is that instead of matching, we formulate the problem as a task of solving an equation. Also, our proposed method only involves a one-step matrix inversion rather than an iterative optimization procedure. The curse of dimensionality may pose some challenges. However, this problem can be addressed by adopting the black-box-shift estimation method, as detailed in the supplementary.

Comment 2 (Questions)

How to select $\alpha$ in terms of the sample size?

Our Response

Thanks for pointing this issue out. We agree that hyperparameter tuning is not trivial. However, based on our Theorem 1, the regularization parameter needs to satisfy the condition that $\lim_{n\to\infty}n\alpha^2=\infty$. Thus, in our experiment, we proposed to choose $\alpha=n^{-1/4}$ as a rule of thumb to meet the condition for practical simplicity. In Section C.2.3 of the supplementary materials, we conducted a sensitivity analysis from which $\alpha=n^{-1/4}$ achieves good performance. We use this rule of thumb for all the numerical studies, and the results are consistently good.

Dear Reviewer 4YDN,

We extend our gratitude for your insightful feedback on our work. Your comments have been invaluable in shaping the refinement of our paper. Here's a summary of the key modifications we made in response to your suggestions:

1. Inclusive comparison: We introduced the L2IWE method proposed by Nguyen et al. (2016) for a comprehensive comparative analysis.
2. Enhanced clarity on novelty and challenge: We redefined the non-trivial target shift problem as a solvable linear equation problem, supported by a robust theoretical foundation.
3. Efficient hyper-parameter tuning: We introduced a simple yet effective rule of thumb method for selecting the regularization parameter, along with a thorough ablation study.

With the end of the rebuttal phase approaching, we kindly invite you to review our responses and the revised paper. Your feedback is crucial to ensuring that our clarifications address the concerns you raised. If our revision and response meet your satisfaction, would you consider raising the score to support our work? We highly value your support and guidance in elevating the quality of our work.

Thank you for your time and consideration.

Authors.

We appreciate your support for our work! Regarding the performance of L2IWE on the SOCR and Szeged Weather datasets, we observed that L2IWE's performance will decrease drastically when the shift becomes more severe. In Section C.3 of our supplementary materials, we presented several examples illustrating the estimated weights for each method. Notably, L2IWE exhibited higher weights concentrated in the middle of the response (i.e., $y$) support, with near-zero weights at the boundaries. However, in these experiments, the oracle weights shifted towards the boundaries. In summary, empirical evidence shows that L2IWE's performance is unsatisfactory under severe shifts. The shifts in the real datasets are more severe; we conjecture that this may contribute to the underperformance of the L2IWE method. Due to the time limit of the discussion period, we do not have enough time to find the causes more in-depth, but we shall investigate this later.

Additionally, we acknowledge the potential impact of ratio density estimation on our algorithm's performance. In Section C.2.3 of the supplementary materials, we conducted a sensitivity study involving two density ratio estimation methods. While we did not observe significant performance changes, we recognize the importance of a more in-depth analysis and discussion on numerical stability. Further exploration in this aspect would contribute to a comprehensive understanding of our algorithm's behavior.

Thanks again for your time and efforts!

Thanks for addressing my comments and implementing the new baseline. The proposed method outperforms the baselines. If possible, it would be useful to discuss why L2IWE fails to beat non-adaptation in the SOCR and Szeged Weather dataset as Table 1 shows negative $\Delta Pred$ MSE. Does any numerical instability occur in the estimation procedure? Then, it would be nice to address why the proposed method does not fail. As other reviewers pointed out, the kernel ratio density estimator could be unstable sometimes. More discussions on numerical stability can be helpful for practitioners to choose between algorithms.

Overall, I think the paper is well-written and the contribution is clear, so I raise my score a bit.

Comment 1

In the experiments, only low-dimensional and small-scale datasets are used for comparison. Thus, the empirical results are not sufficient to support the success of the proposed algorithm. Considering the increasing need of learning with high-dimensional data such as images and the poor performance of kernel density approximation in high-dimensional data, the comparison in high-dimensional data is important.

Our Response

In the revision, we have demonstrated the effectiveness of our proposed method with larger datasets with higher dimensions (e.g., high dimension regression problem on UCI community and crime dataset and knowledge distillation regression problem on MNIST dataset). Due to space limits, we put the additional numerical studies in Sections C.4.1 & 4.2 of the supplementary materials. Furthermore, to handle the density approximation problem with high-dimensional data, we adopted the black-box-shift-estimation method Lipton et al. (2018) proposed in our setting. Please refer to Section C.4 for the details.

Comment 2

Determining the hyperparameter $\alpha$ for different data sets is difficult.

Our Response

Thanks for pointing this issue out. Based on Theorem 1, the regularization parameter must satisfy the condition that $\lim_{n\to\infty}n\alpha^2=\infty$. Thus, in our experiment, we proposed to choose $\alpha=n^{-1/4}$ as a rule of thumb to meet the condition for practical simplicity. In Section C.2.3 of the supplementary materials, we conducted a sensitivity analysis from which $\alpha=n^{-1/4}$ achieves good performance. We use this rule of thumb for all the numerical studies, and the results are good.

Comment 3

In classical ratio estimation approaches, the estimation variance is sometimes large, e.g. $p_s(y)$ goes to zero, making the algorithm unstable. In the continuous label space, it is unclear whether this problem becomes worse. It is suggested to discuss this issue or to introduce an additional variance reduction mechanism.

Our Response

Indeed, it is problematic when $p_s(y)$ goes to zero. This is why we need the assumption that $p_s(y)$ is bounded from below by some constant $\epsilon>0$ (see Assumption 3). Also, we need the support of $p_t(y)$ to be a subset of the support of $p_s(y)$. When $p_s(y)$ goes to zero (while $p_t(y)$ is not zero), it corresponds to the out-of-distribution (OOD) setting. Although this is beyond the scope of our paper, it represents an interesting direction for future research. We discussed this issue as a future research direction in the final discussion.

Comment 4 (Questions)

See the comments in weakness part. As also, it is suggested to test the performance of the proposed algorithm on large-scale and high-dimensional datasets.

In this revision, we applied our method to larger datasets with higher dimensions in additional experiments. More specifically, we include the experiments of high dimension regression problem on the UCI community and crime dataset and knowledge distillation regression problem on the MNIST dataset in the supplementary. The proposed method demonstrated good performance. Due to space limits, the numerical study results can be found in Sections C.4.1 & 4.2 of the supplementary materials.

Dear Reviewer 2aLy,

We express our gratitude for your time and efforts for reviewing our paper and offering valuable suggestions. In our rebuttal, we delved deeper into three key areas:

1. Empirical performance comparison on large-scale and high-dimensional datasets.
2. Elaboration on the regularization parameter selection, along with its ablation study.
3. Further extension of discussions regarding $p_s(y)\rightarrow 0$ and addressing out-of-distribution (OOD) problems.

As we approach the conclusion of the discussion phase, we hope you have had the opportunity to read our responses and the revised paper. In consideration of time constraints, we hope you can inform us of any lingering questions or concerns that require additional attention or clarification. If our revision and response meet your satisfaction, we would greatly appreciate your consideration in raising the score.

Your time and consideration are highly appreciated.

Authors.

Comment 1

Given the importance weight function, is the estimation process sensitive to noise or outliers in the data and thus potentially leading to inaccurate estimates of the importance weight function? Or this limitation is conquered by applying the regularization technique in your paper?

Our Response

Yes, the noise and outlier data will affect the estimation performance and the regularization applied in our method can enhance our ability to robustly adjust to uncertainties in estimating the importance weight function.

Previous studies have illustrated such advantage of regularization in classification tasks; for instance, Azizzadenesheli et al. (2019) introduced a regularized version of the label shift estimator called RLLS, building upon the BBSE method of Lipton et al. (2018), who did not employ regularization. The regularized estimator RLLS demonstrated superior performance to the original non-regularized version BBSE, indicating that the regularization does help.

In the continuous continuous, regularization serves a dual purpose: addressing the ill-posedness problem while providing the advantage of robustness. However, unlike the classification case where regularization is optional (i.e., Lipton et al. 2018), the problem we consider becomes ill-posed due to the 'explosive behavior' in the continuous case without regularization, as detailed in our paper. So, regularization is a must-have in our case.

Comment 2

Estimating a continuous importance weight function is a probabilistic process, together with the regularization process, this may inherent uncertainty in the estimates, affecting the model's reliability in certain cases. How do you solve this issue?

Our Response

The functional estimation, and indeed any estimator, is inherently contingent upon the randomness of the data sample; uncertainty becomes an inherent characteristic of the estimation. In our study, we establish a theoretical bound for the estimator’s convergence rate in Theorem 1, and the error bound measures the estimation uncertainty from the theoretical aspect. Also, in our numerical evaluation, we presented all the estimation performance with confidence intervals and standard deviations, which can capture the estimation uncertainty.

A related aspect involves the uncertainty quantification for the continuous case. For instance, considering the randomness inherent in the estimation, how to construct a prediction set for the response. While this lies beyond the scope of our current paper, it represents an interesting direction for future research. Existing efforts have been made on uncertainty quantification in the classification case (e.g., Podkopaev and Ramdas 2021), but notably, there appears to be a gap in addressing the continuous case. We have added a discussion on future research in the paper.

Comment 3

Regularization can help stabilize the estimation process, but accurate estimation of a continuous importance weight function may still require a substantial amount of data from both source and target domains. What is the sample complexity of your algorithm?

Our Response

We agree that the sample sizes play an important role in the weight estimation. Theorem 1 provides a stochastic order bound for our proposed functional estimation, in which the sample size will directly affect the convergence rate. In Equation (11), the bound is a function of $n$ (sample size of the source data) and $m$ (sample size of the target data) given other quantities (e.g., the tuning parameter $\alpha$ and the bandwidth $h$).

References

Lipton, Z., Wang, Y. & Smola, A.. (2018). Detecting and Correcting for Label Shift with Black Box Predictors. ICML.

Azizzadenesheli, K., Liu, A., Yang, F., & Anandkumar, A. (2019). Regularized learning for domain adaptation under label shifts. ICLR.

Podkopaev, A., & Ramdas, A. (2021). Distribution-free uncertainty quantification for classification under label shift. UAI.

Dear Reviewer D6fP,

We would like to thank you for your thorough review and insightful suggestions. In response to your feedback, we have incorporated additional discussion points in our rebuttal:

1. Explanation of how regularization contributes to estimation robustness.
2. Detailed insights into the measurement and quantification of estimation uncertainty.
3. Clarification of the theoretical sample complexity of our proposed method.

As we near the conclusion of the rebuttal phase, we kindly invite you to take a moment to review our responses. We are eager to know if our efforts have addressed and clarified the concerns you raised. If our revision and response align with your satisfaction, we would sincerely appreciate your consideration in raising the score.

Your time and consideration are greatly appreciated.

Thank you,

Authors.