Fine-Grained Dynamic Framework for Bias-Variance Joint Optimization on Data Missing Not at Random

Response for Weaknesses A1. For the discussion of figures, the descriptions of theoretical results and their proofs, and experimants, according to comments from reviewers, we will add more explanation to avoid any ambiguity and improve the readability.

A2. Sorry for the unprecise descriptions. We will modify the descriptions in lines 99-100 to the following descriptions.

''..., then IPS and DR estimators are unbiased. For a new dataset, we cannot know in advance the range of the propensities in this dataset. Therefore, a new dataset may introduce extremely small propensities to lead to unbounded variances of IPS and DR, which will disrupt the stability of estimators, especially for larger datasets. It is unacceptable for real industrial scenarios.''

A3. Considering (1), for all $(u,i)$, we have $f(o _{u,i},\hat{p} _{u,i})e _{u,i}+g(o _{u,i},\hat{p} _{u,i})\hat{e} _{u,i}\ge0$ and $h(o _{u,i},\hat{p} _{u,i})\ge0$. Let $f(o _{u,i},\hat{p} _{u,i})e _{u,i}+g(o _{u,i},\hat{p} _{u,i})\hat{e} _{u,i}$ be denoted as $r(o _{u,i},\hat{p} _{u,i},e _{u,i},\hat{e} _{u,i})$. As given in the proof of Theorem 3.1, we have

$ \text{Cov}(L _\text{Est}, L _\text{Reg})\=\frac{1}{\vert\mathcal{D}\vert^2}\sum _{j=1}^{\vert{D}\vert}\sum _{k=1}^{\vert{D}\vert}\mathbb{E} _O[h(o _j,\hat{p} _j)r(o _k,\hat{p} _k,e _k,\hat{e} _k)]. $

Since $\mathbb{E} _O[h(o _j,\hat{p} _j)r(o _k,\hat{p} _k,e _k,\hat{e} _k)]\ge0$, we obtain $\text{Cov}(L _{Est}, L _{Reg})=\mathbb{E} _O(L _\text{Est}L _\text{Reg})\ge0$.

We will add the above formula to proof of Theorem 3.1.

A4. The baseline approaches and the dynamic estimators are conducted on the KuaiRec dataset, where the experiment setting is same as the setting of other datasets. The corresponding results are given as follows:

In the final version, we will merge this table into Table 2.

A5. Thanks for your careful reading. We will correct the language issues and further improve the presentation in the final version.

Response for Questions:

A1. Sorry for the incorrect formula in line 757 and thanks very much for your careful reading. $'='$ should be $'\ge'$ in this formula.

Denote $f(o _{u,i},\hat{p} _{u,i})e _{u,i}+g(o _{u,i},\hat{p} _{u,i})\hat{e} _{u,i}$ as $r(o _{u,i},\hat{p} _{u,i},e _{u,i},\hat{e} _{u,i})$. Then, the term in line 757 satisfies

$ \frac{1}{\vert\mathcal{D}\vert^2}&\mathbb{E} _O\Bigg[\bigg(\sum _{(u,i)\in\mathcal{D}}r(o _{u,i},\hat{p} _{u,i},e _{u,i},\hat{e} _{u,i})\bigg)^2\Bigg]\\\\ =&\frac{1}{\vert\mathcal{D}\vert^2}\sum _{j=1}^{\vert\mathcal{D}\vert}\sum _{k=1}^{\vert\mathcal{D}\vert}\mathbb{E} _O\Big[r(o _j,\hat{p} _j,e _j,\hat{e} _j)r(o _k,\hat{p} _k,e _k,\hat{e} _k)\Big]\\\\ \ge&\frac{1}{\vert\mathcal{D}\vert^2}\sum _{j=1}^{\vert\mathcal{D}\vert}\mathbb{E} _O\big[r^2(o _j,\hat{p} _j,e _j,\hat{e} _j)\big]\quad\quad\quad\text{becase of } r(o _{u,i},\hat{p} _{u,i},e _{u,i},\hat{e} _{u,i})\ge0 $

Therefore, $\mathbb{V}_O[L]$ satisfies

$ \mathbb{V} _O[L]\ge&\frac{1}{\vert\mathcal{D}\vert^2}\sum _{(u,i)\in\mathcal{D}}\frac{1-p _{u,i}}{p _{u,i}}\big[e _{u,i}-g(0,\hat{p} _{u,i})\hat{e} _{u,i}\big]^2. $

Therefore, for any $e _{u,i}-g(0,\hat{p} _{u,i})\hat{e} _{u,i}\ne0$, $\lim _{p _{u,i}\to0}\mathbb{V} _O[L]=\infty$.

We will use the above proof to replace the formulas in lines 757 and 758.

A2. In the optimization problem (4), $E _B(\cdot)$ and $E _V(\cdot)$ are the measure metrics of bias and variance, respectively. As mentioned in Bias-Variance Quantitative Joint Optimization, (4) can be simplified as $w _1h _B^{Est}(\alpha _{u,i})+w _2h _V^{Est}(\alpha _{u,i})$. We can obtain the analytical solution of the optimal parameter $\alpha _{u,i}^{opt}$ given in Theorem 3.5.

A3. Considering the principle Isotonic Propensity, the function $f(\hat{p} _{u,i})$ is a probability mapping. Therefore, $f(\hat{p} _{u,i})$ should be a monotonically increasing function with $f(0)=0$ and $f(1)=1$. $f(\hat{p} _{u,i})>\hat{p} _{u,i}$ guarantees $h _B^{Est}\ge0$ and $h _V^{Est}\ge0$. With this operation, (3) can be simplified as $\min _{\alpha _{u,i}}\{w _1h _B(\alpha _{u,i})+w _2h _V(\alpha _{u,i})\}$. Then, we can obtain the analytical solution of $\alpha _{u,i}^{opt}$ without increasing the computational complexity. Next, $\lim\limits _{\hat{p} _{u,i}\to0}\frac{\hat{p} _{u,i}}{f(\hat{p} _{u,i})}=C,\forall\alpha _{u,i}\in[0,1]$ ensures

$ \lim _{\hat{p} _{u,i}\to0}h^{\text{Est}} _B(\hat{p} _{u,i},p _{u,i},\alpha _{u,i})=1-Cf^{1-\alpha _{u,i}}(0) $

and

$ \lim _{\hat{p} _{u,i}\to0,\alpha _{u,i}\leq0.5}h^{\text{Est}} _V(\hat{p} _{u,i},p _{u,i},\alpha _{u,i})=\lim _{\hat{p} _{u,i}\to0,\alpha _{u,i}\leq0.5}C(1-p _{u,i})f^{1-2\alpha _{u,i}}(\hat{p} _{u,i})=0 $

In final version, we will add the above discussion to explain these two function design principles.

A1. Thanks for your valuable comments. In the final version, we will add more metrics. Besides, the baseline approaches and the dynamic estimators are conducted on the KuaiRec dataset to verify the performance of the developed dynamic fine-grained framework. The experimental results will be merged into Table 1.

A2. Based on $\alpha _{u,i}\in[0,1]$, when $w _2/w _1$ satisfies $\frac{w _2}{w _1}\leq\frac{f(p _{u,i})}{2(1-p _{u,i})}$, $\alpha _{u,i}^{opt}=1$; when $\frac{f(p _{u,i})}{2(1-p _{u,i})}\leq\frac{w _2}{w _1}\leq\frac{1}{2(1-p _{u,i})}$, then $\alpha _{u,i}^{opt}=\frac{\ln\Big(\frac{2w _2}{w _1}(1-p _{u,i})\Big)}{\ln(f(p _{u,i}))}$; when $\frac{w _2}{w _1}\ge\frac{1}{2(1-p _{u,i})}$, then $\alpha _{u,i}^{opt}=0$. Therefore, the daynamic estimator can be rewritten as

$ L _{D-IPS}=&\frac{1}{\vert\mathcal{D}\vert}\sum _{(u,i)\in\mathcal{D}}\mathbb{1} _{\bigg[\frac{w _2}{w _1}\leq\frac{f(p _{u,i})}{2(1-p _{u,i})}\bigg]}\frac{o _{u,i}}{f(\hat{p} _{u,i})}e _{u,i}+\mathbb{1} _{\bigg[\frac{f(p _{u,i})}{2(1-p _{u,i})}\leq\frac{w _2}{w _1}\leq\frac{1}{2(1-p _{u,i})}\bigg]}\frac{o _{u,i}}{f^{\alpha^{opt} _{u,i}}(\hat{p} _{u,i})}e _{u,i} +\mathbb{1} _{\bigg[\frac{w _2}{w _1}\ge\frac{1}{2(1-p _{u,i})}\bigg]}o _{u,i}e _{u,i},\\\\ L _{D-DR}=&\frac{1}{\vert\mathcal{D}\vert}\sum _{(u,i)\in\mathcal{D}}\mathbb{1} _{\bigg[\frac{w _2}{w _1}\leq\frac{f(p _{u,i})}{2(1-p _{u,i})}\bigg]}\bigg(\hat{e} _{u,i}+\frac{o _{u,i}}{f(\hat{p} _{u,i})}\delta _{u,i}\bigg) +\mathbb{1} _{\bigg[\frac{f(p _{u,i})}{2(1-p _{u,i})}\leq\frac{w _2}{w _1}\leq\frac{1}{2(1-p _{u,i})}\bigg]}\bigg(\hat{e} _{u,i}+\frac{o _{u,i}}{f^{\alpha^{opt} _{u,i}}(\hat{p} _{u,i})}\delta _{u,i}\bigg) +\mathbb{1} _{\bigg[\frac{w _2}{w _1}\ge\frac{1}{2(1-p _{u,i})}\bigg]}\bigg(\hat{e} _{u,i}+o _{u,i}\delta _{u,i}\bigg). $

Therefore, the proposed framework dynamically selects the most appropriate estimator for each user-item pair.

$f(\hat{p} _{u,i})>\hat{p} _{u,i}$ guarantees $h _B^{Est}\ge0$ and $h _V^{Est}\ge0$. With this operation, the joint optimization problem (3) can be simplified as $\min _{\alpha _{u,i}}\{w _1h _B(\alpha _{u,i})+w _2h _V(\alpha _{u,i})\}$. Then, we can obtain the analytical solution of $\alpha _{u,i}^{opt}$ without increasing the computational complexity.

A3. In (5), we have given the analytical solution of the optimal parameter. If an engineer have an IPS-based or DR-based estimators for large-scale datasets, then she/he can modify the IPS-based or DR-based estimators into a D-IPS-based or D-DR-based estimator by the following core code.

w1 = 1
w2 = 0.5
star = 2 * w2 * (1 - propensity) / w1
alpha = np.log(star) / np.log(propensity) 
lower_bound = np.zeros(np.size(propensity)) 
upper_bound = np.ones(np.size(propensity))
alpha = np.where(alpha > upper_bound, upper_bound, alpha)
alpha = np.where(alpha < lower_bound, lower_bound, alpha)
propensity = np.power(propensity, alpha) 

A4. In the proposed dynamic framework, the weight ratios and the function forms jointly determined the performance of estimators. As mentioned in A2, when $\frac{w _2}{w _1}=0$ and $f(\hat{p} _{u,i})=\hat{p} _{u,i}$, then the optimal parameter $\alpha _{u,i}^{opt}$ is set as 1, and the D-IPS and D-DR methods are equivalent to IPS and DR methods, which are unbiased approaches. According to Fig. 2, we can obeserved the performance of dynamic estimators under different weight ratio, which includes the case of $\frac{w _2}{w _1}=0$. On the other hand, under the identical weight ratio $\frac{w _2}{w _1}=0.1$, effects of different functions on performances of dynamic estimators are shown in Table 3. From $Gain _{AUC}$ and $Gain _{N}$ in Table 3, we can observe the individual contribution of different function forms.

A5. The definition of the variance reveals the preformance stability of the estimator in whole propensity space. For a new dataset, we cannot know in advance the range of the propensities in this dataset. An unknown biased dataset may disrupt the stability of the prediction model, which brings significant risks to practical applications. Therefore, it is necessary to consider the boundedness of variance and generalization to ensure the stability and performance of prediction models. In real-world applications, the theoretical bounds for variance and generalization indrectly reflect the stability and generalization performance of estimators, respectively.

A6. We find that the function $f(\hat{p} _{u,i})$ is not unique. According to the functon design principles, we can obtain a family of functions. According to Fig. (c), for a fixed propensity, different functions correspond to different minimum objective values. For different datasets with different propensity distributions, we can determine the function form according to the minimum objective values of functions under the propensity distribution. Actually, there is further optimization potential for the selection of the function form. We can incorporate the selection of function forms into the bias-variance joint optimization problem (4). The corresponding optimization problem can be rewritten as ${Objective}^{opt}=\min _{\alpha _{u,i}\in[0,1], f\in\mathcal{F}}\Big[w _1E _B(h^{Est} _B(f^{\alpha _{u,i}}(\hat{p} _{u,i})))+w _2E _V(h^{Est} _V(f^{\alpha _{u,i}}(\hat{p} _{u,i})))\Big],{s.t.} 0\leq\alpha _{u,i}\leq1,$ where $\mathcal{F}$ is the candidate function set. This problem is one of our future works.

A7. $w _1$ and $w _2$ in Eq. (3) are the hyper-parameter. $\alpha _{u,i}^{opt}$ only depends on $w _2/w _1$. $w _2/w _1$ will leads to the Pareto frontier of the estimation performance. As of now, optimization problem for weight ratios is still an open question because of unknown properties of its Pareto frontier. In practice, for a new dataset, the tuning range of $w _2/w _1$ is usually set to be between 0 and 1. Then several $w _2/w _1$ are used to test the estimation performance and the optimal weight ratio is selected as the final parameter.

A1. Considering (1), for all $(u,i)$ pairs, we have $f(o _{u,i},\hat{p} _{u,i})e _{u,i}+g(o _{u,i},\hat{p} _{u,i})\hat{e} _{u,i}\ge0$ and $h(o _{u,i},\hat{p} _{u,i})$. To facilitate representation, $f(o _{u,i},\hat{p} _{u,i})e _{u,i}+g(o _{u,i},\hat{p} _{u,i})\hat{e} _{u,i}$ is denoted as $r(o _{u,i},\hat{p} _{u,i},e _{u,i},\hat{e} _{u,i})$. In the proof of Theorem 3.1, $\text{Cov}(L _\text{Est}, L _\text{Reg})$ satisfies

$ \text{Cov}(L _\text{Est}, L _\text{Reg})=&\frac{1}{\vert\mathcal{D}\vert^2}\mathbb{E} _O\Bigg(\sum _{(u,i)\in\mathcal{D}}\bigg[h(o _{u,i},\hat{p} _{u,i})\sum _{(u,i)\in\mathcal{D}}r(o _{u,i},\hat{p} _{u,i},e _{u,i},\hat{e} _{u,i})\bigg]\Bigg)\nonumber\\\\=&\frac{1}{\vert\mathcal{D}\vert^2}\mathbb{E} _O\Bigg(\sum _{j=1}^{\vert{D}\vert}\sum _{k=1}^{\vert{D}\vert}h(o _j,\hat{p} _j)r(o _k,\hat{p} _k,e _k,\hat{e} _k)\Bigg)\nonumber\\\\=&\frac{1}{\vert\mathcal{D}\vert^2}\sum _{j=1}^{\vert{D}\vert}\sum _{k=1}^{\vert{D}\vert}\mathbb{E} _O[h(o _j,\hat{p} _j)r(o _k,\hat{p} _k,e _k,\hat{e} _k)]. $

Since $\mathbb{E} _O[h(o _j,\hat{p} _j)r(o _k,\hat{p} _k,e _k,\hat{e} _k)]\ge0$, we obtain $\text{Cov}(L _{Est}, L _{Reg})\ge0$.

We will revise this sentense in line 124 as follows, and add the above formula (2) to proof of Theorem 3.1.

"... respectively. For all $(u,i)$ pairs, they satisfy $f(o _{u,i},\hat{p} _{u,i})e _{u,i}+g(o _{u,i},\hat{p} _{u,i})\hat{e} _{u,i}\ge0$ and $h(o _{u,i},\hat{p} _{u,i})$. $\lambda>0$ is a ..."

Since Corollary 3.2 is the contrapositive of Theorem 3.1, the proof of Theorem 3.1 has been given. Therefore, Corollary 3.2 also holds. We further provide the proof of Corollary 3.2 as follows

proof of Corollary 3.2

We use the method of proof by contradiction. Assume that when $\mathbb{V} _O[L _\text{Est+Reg}]\leq\mathbb{V} _O[L _\text{Est}]$, $L _\text{Est+Reg}$ is unbiased.

According to $\mathbb{V} _O[L _\text{Est+Reg}]\leq\mathbb{V} _O[L _\text{Est}]$, we have

$ \mathbb{V} _O[L _\text{Est+Reg}]=&\mathbb{V} _O[L _\text{Est}]+2\lambda\text{Cov}(L _\text{Est}, L _\text{Reg})+\lambda^2\mathbb{V} _O[L _\text{Reg}]\nonumber\\\\ \leq&\mathbb{V} _O[L _\text{Est}], $

which implies that $2\lambda\text{Cov}(L _\text{Est}, L _\text{Reg})+\lambda^2\mathbb{V} _O[L _\text{Reg}]\leq0$. Therefore, the parameter $\lambda$ needs to satisfy $0\leq\lambda\leq-\frac{2\text{Cov}(L _\text{Est}, L _\text{Reg})}{\mathbb{V} _O[L _\text{Reg}]}$, which implies $\text{Cov}(L _\text{Est}, L _\text{Reg})\leq0$.

As shown in A1, when $L _\text{Est+Reg}$ is unbiased, we have $\text{Cov}(L _{Est}, L _{Reg})\ge0$, which contradicts the condition $\text{Cov}(L _\text{Est}, L _\text{Reg})\leq0$. Therefore, Corollary 3.2 holds.

A2. We can obtain a family of functions. According to $w _1h _B^{Est}(\alpha,p _{u,i})+w _2h _V^{Est}(\alpha,p _{u,i})$ and Fig. (c), we can observe that different $f^{\alpha _{u,i}}(\hat{p} _{u,i})$ will lead to different objective function surface. For different datasets with different propensity distributions, we can determine the function form according to the minimum objective values of functions. Actually, there is further optimization potential for the selection of the function form. We can incorporate the selection of function forms into the bias-variance joint optimization problem (4), which is formulated as ${Objective}^{opt}=\min _{\alpha _{u,i}\in[0,1], f\in\mathcal{F}}\Big[w _1E _B(h^{Est} _B(f^{\alpha _{u,i}}(\hat{p} _{u,i})))+w _2E _V(h^{Est} _V(f^{\alpha _{u,i}}(\hat{p} _{u,i})))\Big],{s.t.} 0\leq\alpha _{u,i}\leq1$ This problem is one of our future works.

A3. This work focus on revealing the general rules of regularization techniques and unbiased estimators and developed a fine-grained dynamic framework to jointly optimize bias and variance. Therefore, almost all learning algorithms can be transformed into the fine-grained dynamic framework dynamic framework, such as joint learning of DR, double learning of MRDR, cycle learning of SDR, collaborative learning of TDR. In the final version, we will add the following core code of fine-grained dynamic estimators to Appendix E to explain how to transform a propensity-based debiased approach to the fine-grained dynamic estimator.

Core Code of Fine-Grained Dynamic Estimators

w1 = 1
w2 = 0.5
rate = w2 / w1
star = 2 * rate * (1 - propensity)
# corresponds to function 1 given in Table 1
alpha = np.log(star) / np.log(propensity) 
# alpha = np.log(star) / np.log(np.sin(propensity)/np.sin(1)) 
# alpha = np.log(star) / np.log(np.log(propensity+1)/np.log(2)) 
# alpha = np.log(star) / np.log(np.tanh(propensity)/np.tanh(1)) 
lower_bound = np.zeros(np.size(propensity))
upper_bound = np.ones(np.size(propensity))
alpha = np.where(alpha > upper_bound, upper_bound, alpha)
alpha = np.where(alpha < lower_bound, lower_bound, alpha)
propensity = np.power(propensity, alpha) 

A4. In the final version, we will add the following descriptions in line 313 to improve Related Work.

“...forth. To improve the robustness of estimators, a multiple robust estimator is developed in [3] by taking the advantage of multiple candidate imputation and propensity models, which is unbiased when any of the imputation or propensity models, or a linear combination of these models is accurate. From a novel function balancing perspective, Li et al. propose to approximate the balancing functions in reproducing kernel Hilbert space and design adaptively kernel balancing IPS/DR learning algorithms [1]. Moreover, aimed at limitations of miscalibrated imputation and propensity models, Kweon and Yu [2] propose a doubly calibrated estimator and a tri-level joint learning framework to simultaneously optimize calibration experts alongside prediction and imputation models. For the...”

A5. Thanks for your careful reading. We will correct typos, language issues, and unclear desciption, and further improve the presentation in the final version.

Response for Weaknesses:

• In the final version, we will add the following analysis of Figure 1 to lines 212 and 220, respectively.

''Line 212:.... The curves of objective functions under different designed functions $f(\cdot)$ are given in Fig. 1(c). It can be oberved that for a fixed propensity, there exists an $\alpha$ such that the objective function attains the minimum value. Besides, ...''

''Line 220:.... Under different designed function $f(\cdot)$, the schematic diagram of optimal objective values corresponding to the optimal parameter $\alpha_{u,i}^{opt}$ is shown in Figure 1(d). Next, ...''

• The constant $C$ is not related to $\alpha_{u,i}$. In line 189, what we want to convey is that for all $\alpha_{u,i}\in[0,1]$, the constant $C$ is the same.

• Yes, $w_1$ and $w_2$ in Eq. (3) are the hyper-parameter. As shown in (5), $\alpha_{u,i}^{opt}$ depends on the weight ratio $w_2/w_1$. The weight ratio space will leads to the Pareto frontier of prediction model performance. As disussed in Conclusions, various properties of the Pareto frontier and optimization methods for weight ratios are still an open question, and it is one of our future works. In practice, the tuning range of $w_2/w_1$ is usually set to be between 0 and 1.

• In Theorem 3.5, when $\frac{w_2}{w_1}=0$ and $f(\hat p_{u,i})=\hat p_{u,i}$, then the optimal parameter $\alpha_{u,i}^{opt}$ is set as 1, and D-IPS and D-DR are equivalent to IPS and DR, which are unbiased. Therefore, the unbiased conditions are $\frac{w_2}{w_1}=0$, $f(\hat p_{u,i})=\hat p_{u,i}$ when propensities or imputation errors are accurate. When $\alpha$ is optimal, according to the Lemma D.1, the bias of the proposed estimator can be obtained. When $\frac{w_2}{w_1}\leq\frac{f(p_{u,i})}{2(1-p_{u,i})}$, then $\alpha_{u,i}^{opt}=1$; when $\frac{f(p_{u,i})}{2(1-p_{u,i})}\leq\frac{w_2}{w_1}\leq\frac{1}{2(1-p_{u,i})}$, then $\alpha_{u,i}^{opt}=\frac{\ln\Big(\frac{2w_2}{w_1}(1-p_{u,i})\Big)}{\ln(f(p_{u,i}))}$; when $\frac{w_2}{w_1}\ge\frac{1}{2(1-p_{u,i})}$, then $\alpha_{u,i}^{opt}=0$. Therefore, the dynamic estimators are applicable to different propensity distributions by adjusting the weight ratio $w _2/w _1$. Different weight ratios will lead to different distributions of instances on different estimators. Therefore, even if there is no propensity score that will not tend to 0, the developed fine-grained dynamic framework still can achieve quantitative optimization of bias and variance.

• Equation (10) should be

$ \mathbb{V} _O[L _\text{Est+Reg}] =\mathbb{E} _O[L^2 _\text{Est}+\lambda^2{L}^2 _\text{Reg}+2\lambda{L} _\text{Est}L _\text{Reg}]-\mathbb{E}^2 _O[L _\text{Est}+\lambda{L} _\text{Reg}]. $

Response for Questions:

A1. Considering (1), for all $(u,i)$ pairs, we have $f(o _{u,i},\hat{p} _{u,i})e _{u,i}+g(o _{u,i},\hat{p} _{u,i})\hat{e} _{u,i}\ge0$ and $h(o _{u,i},\hat{p} _{u,i})\ge0$. To facilitate representation, $f(o _{u,i},\hat{p} _{u,i})e _{u,i}+g(o _{u,i},\hat{p} _{u,i})\hat{e} _{u,i}$ is denoted as $r(o _{u,i},\hat{p} _{u,i},e _{u,i},\hat{e} _{u,i})$. As given in the proof of Theorem 3.1, $\text{Cov}(L _\text{Est}, L _\text{Reg})$ satisfies

$ \text{Cov}(L _\text{Est}, L _\text{Reg})=&\frac{1}{\vert\mathcal{D}\vert^2}\mathbb{E} _O\Bigg(\bigg[\sum _{(u,i)\in\mathcal{D}}r(o _{u,i},\hat{p} _{u,i},e _{u,i},\hat{e} _{u,i})\bigg]\bigg[\sum _{(u,i)\in\mathcal{D}}h(o _{u,i},\hat{p} _{u,i})\bigg]\Bigg)\nonumber\\\\=&\frac{1}{\vert\mathcal{D}\vert^2}\mathbb{E} _O\Bigg(\sum _{(u,i)\in\mathcal{D}}\bigg[h(o _{u,i},\hat{p} _{u,i})\sum _{(u,i)\in\mathcal{D}}r(o _{u,i},\hat{p} _{u,i},e _{u,i},\hat{e} _{u,i})\bigg]\Bigg)\nonumber\\\\=&\frac{1}{\vert\mathcal{D}\vert^2}\mathbb{E} _O\Bigg(\sum _{j=1}^{\vert{D}\vert}\sum _{k=1}^{\vert{D}\vert}h(o _j,\hat{p} _j)r(o _k,\hat{p} _k,e _k,\hat{e} _k)\Bigg)\nonumber\\\\=&\frac{1}{\vert\mathcal{D}\vert^2}\sum _{j=1}^{\vert{D}\vert}\sum _{k=1}^{\vert{D}\vert}\mathbb{E} _O[h(o _j,\hat{p} _j)r(o _k,\hat{p} _k,e _k,\hat{e} _k)]. $

Since $\mathbb{E} _O[h(o _j,\hat{p} _j)r(o _k,\hat{p} _k,e _k,\hat{e} _k)]\ge0$, we obtain $\text{Cov}(L _{Est}, L _{Reg})=\mathbb{E} _O(L _\text{Est}L _\text{Reg})\ge0$.

We will revise this sentense in line 124 as follows, and add the above formula (1) to proof of Theorem 3.1.

... respectively. For all $(u,i)$ pairs, they satisfy $f(o _{u,i},\hat{p} _{u,i})e _{u,i}+g(o _{u,i},\hat{p} _{u,i})\hat{e} _{u,i}\ge0$ and $h(o _{u,i},\hat{p} _{u,i})\ge0$. $\lambda>0$ is a ...

A2. 'the one' in Theorem 3.1 and Corollary 3.2 is a pronoun that refers to the variance. What we want to convey in Theorem 3.1 and Corollary 3.2 are that

"Theorem 3.1 Let $L _\text{Est+Reg}$ be defined in (1) and the estimator $L _\text{Est}$ be unbiased. If $L _\text{Est+Reg}$ is unbiased, then the variance of $L _\text{Est+Reg}$ is greater than the variance of the original estimator $L _\text{Est}$."

"Corollary 3.2 If the variance of $L _\text{Est+Reg}$ is less than the variance of the original estimator $L _\text{Est}$, then $L _\text{Est+Reg}$ is not unbiased."

In the final version, we will modify Theorem 3.1 and Corollary 3.2 to the descriptions above to avoid any ambiguity.

A3. According to the proof of Theorem 3.3, we can obtain

$ \mathbb{V} _O[L _\text{Est+Reg}]\ge\mathbb{E} _O[L^2 _\text{Est}]+\lambda^2\mathbb{E} _O[{L}^2 _\text{Reg}]+2\lambda\mathbb{E} _O[L _\text{Est}L _\text{Reg}]-\bar{B}^2.\nonumber $

Based on A1., we can obtain $\lambda^2\mathbb{E} _O[{L}^2 _\text{Reg}]\ge0$ and $\text{Cov}(L _{Est}, L _{Reg})\ge0$. Since $\mathbb{E} _O[L^2 _\text{Est}]$ tends to infinity, $\mathbb{E} _O[L^2 _\text{Est}]+\lambda^2\mathbb{E} _O[{L}^2 _\text{Reg}]+2\lambda\mathbb{E} _O[L _\text{Est}L _\text{Reg}]$ also tends infinity. For all estimators, ${L} _\text{Est}\ge0$ and ${L} _\text{Reg}\ge0$, ${L} _\text{Est}=-{L} _\text{Reg}$ is true only when ${L} _\text{Est}={L} _\text{Reg}=0$.