First, thanks for careful reading our submission. Please let us reply to your questions one by one.

[Q1 The fundamental difference between (conditional) CS divergence and $I(y;t)$ with Gaussian variational assumption]

In fact, “CS+Gaussian” and “KL+Gaussian” are fundamentally different. Please refer to our reply to your [Q2].

We now reply to your second question "why the CS-IB outperforms others also for $r=0$". In fact, when the compression ratio $r=0$, all competing methods do not have the compression term, our CS-IB optimizes $D_{CS}(p(y|x);q_\theta(\hat{y}|x))$ which reduces to MSE when kernel size $\sigma=0$; whereas other methods optimize MSE directly. For our method, we never select kernel size $\sigma=0$. This is because by our analysis in Appendix B.4 in the revised submission, $\sigma=0$ significantly increases the estimation variance (which make the results highly unstable).

There are two reasons why we are better: 1) we do not make parametric assumption on $q_\theta(\hat{y}|x)$. According to our analysis in Appendix B.2, a parametric assumption like Gaussian might not be the optimal choice, and the estimation bias is somewhat hard to control. In fact, our performance is better further justifying that there is a large bias introduced by the Gaussian assumption. Note that, for variational KL divergence estimator, the bias is exactly $D_{KL}(p(y|x);q_\theta(\hat{y}|x))$ or $D_{KL}(p(y|t);q_\theta(\hat{y}|t))$ (please refer to Eq.(96) in Appendix B.2).

Another reason why $D_{KL}(p(y|x);q_\theta(\hat{y}|x))$ performs better than MSE when $\sigma\neq0$ is because, CS divergence encourages that if two points $x_i$ and $x_j$ are sufficiently close, their predictions ($\hat{y}_i$ and $\hat{y}_j$) will be close as well, which essentially enhances the numerical stability of the trained model. When $\sigma=0$, there is no such regularization.

[Q2 The inductive bias between KL+Gaussian and CS+Gaussian]

In fact, “KL+Gaussian” and “CS+Gaussian” are fundamentally different: the former is a parametric estimator that assumes the underlying distribution is Gaussian; whereas the latter is a non-parametric estimator that does not make any assumption on underlying distribution. A Gaussian kernel in kernel density estimator (KDE) does not mean we assume the underlying data distribution is Gaussian.

When the underlying distribution is not Gaussian, “CS+Gaussian” still works, whereas “KL+Gaussian” suffers from poor performance. We added a toy example in Appendix C.5 to illustrate this point.

During the rebuttal, we also systematically analyze the bias and variance of our non-parametric estimators and show that they are consistent with a rate of convergence of $N\sqrt{\sigma}$, in which $N$ is the number of samples, $\sigma$ is the kernel width. Please refer to our discussion in Appendix B.4.

Moreover, we also analyze the bias of previous variational KL divergence-based methods. Please refer to our discussion in Appendix B.2. Our argument is that previous literature either has obvious much more bias than ours or similar to us.

To summarize, for our estimator, the asymptotic property is guaranteed. But the use of Gaussian kernel in KDE may introduce some inductive bias, since other kernels may perform better (e.g., smaller variance) or more robust, such as Epanechnikov kernel. But for “KL+Gaussian”, the Gaussian assumption alone for the underlying distribution makes the estimation is highly likely to be biased.

[Q3 substantial mathematical advantages by CS divergence]

Thanks for the good question. First, we agree with you that in a deterministic network, the true CS-QMI between $x$ and $t$ may also become infinity. To check why, let us consider the CS-QMI between two Gaussian signals $x$ and $t$ with correlation coefficient $\rho$, in which $I_{CS}(x;t)$ is given by: $I_{CS}(x;t) = \frac{1}{2}\log\left( \frac{4-\rho^2}{4\sqrt{1-\rho^2}} \right)$. Suppose $x$ passes through an identity mapping, that is $t=x$, and $\rho=1$, we have $I_{CS}(x;t)\rightarrow\infty$.

However, we would like to make a clarification that the CS divergence is not a special case of the well-known Renyi’s divergence developed by Renyi with the following expression: D_\alpha$$(p;q) = \frac{1}{\alpha-1}\log\int p^\alpha q^{1-\alpha}du=\frac{1}{\alpha-1}\log\int p\left(\frac{q}{p}\right)^{1-\alpha}du=\frac{1}{\alpha-1}\log\mathbb{E}_p\left(\frac{q}{p}\right)^{1-\alpha}

When $\alpha=2$, we obtain a new divergence $-\log(1-\chi^2(p;q))$ that is proportional to the $\chi^2$ divergence, rather than our CS divergence mentioned in the paper.

To our knowledge, there are two big ideas to define divergence: the integral probability metrics (IRM) aim to measure $p-q$, whereas the Renyi’s divergences with order $\alpha$ aim to model the ratio of $\frac{p}{q}$. The general idea of CS divergence is totally different: it neither horizontally computes $p-q$ nor vertically evaluates $\frac{p}{q}$. Rather, it quantifies the “distance” of two distributions by quantify the tightness (or gap) of an inequality associated with integral of densities (i.e., $\int pq$). That is, CS divergence does not belong to either IRM or the traditional Renyi’s divergence family. Please refer to our discussion on Appendix B.1 for more details.

This uniqueness makes CS divergence has many desirable properties compared to KL. One of them is that it relaxes the constraint on the support. Note that, $D_{KL}(p;q)$ has finite values if and only if support of $p$ is a subset of support of $q$. This condition is very strong, which also explains why we use KL divergence combined with non-parametric estimator, the model is hard to train (please also refer to our discussion in Appendix B.3, especially the last paragraphs marked with blue).

Finally, we would also like to mention two important and special properties of CS divergence over KL divergence (or Renyi divergence), although these two properties are not explicated used in our method, but it would be beneficial and helpful to the community. (1) the CS divergence has closed-form expression for mixture-of-Gaussians (MoG); (2) CS divergence is invariant to scaling. That is, given $\lambda_1,\lambda_2>0$, we always have $D_{\text{CS}} (p;q) = D_{\text{CS}} (\lambda_1 p;\lambda_2 q)$. Please also refer to our Appendix D for details.

[Q4 Comparison of methods in Table 2 and Fig. 2]

We would like to clarify that in Fig.2, the $y$-axis that quantities $I(y;t)$ is indeed measured with a common currency (to provide a fair comparison). Specifically, we follow (Kolchinsky et al., 2019b), and measure the true value of $I(y;t)$ by $1/2 \log⁡(\text{var}(y)/MSE)$. Hence, $I(y;t)$ in Fig.2 is not measured with respective estimators of different methods.

Going back to the Table 2, yes, the compression ratio $r$ is indeed computed with respective estimators. However, we would like to clarify here that $r$ is a relative quantity, rather than an absolute quantity. Specifically, $r$ is defined as: $r = 1- \frac{a}{b}$, in which $a=I(x;t)$ when $\beta=\beta^*$ and $b=I(x;t)$ when $\beta=0$ (i.e., no compression).

Hence, for each method, the physical meaning of $r$ is how much information has been compressed (when $\beta=\beta^*$) relatively to no compression (i.e., $\beta=0$).

We use $r$ rather than $I(x;t)$ in both Table 2 and Fig 2. This is just because different IB approaches differ in the way to (define) or estimate $I(x;t)$. Most literatures use Shannon’s mutual information (MI), our CS-QMI estimates a generalized MI. By contrast, HSIC-bottleneck uses HSIC, which does not have physical meaning of mutual information in terms of bits or nats. Hence, if we use $I(x;t)$, different approaches have different units, which may result in an unfair comparison.

Finally, we would like to emphasize that, for each method, we additionally report the best performance achieved when $r\neq0$ (i.e., $β\neq0$) in Table 7 in Appendix C.4. This result suggests that, under an optimal compression ratio (in which the $r$ or $\beta$ is different for different methods), our CS-IB is significantly and consistently better.

First, thanks for your positive comments. Please let us reply to your questions one by one.

[Q1 A generalization of Cauchy-Schwarz divergence]

Thanks for the good question. Yes, you are right, the general idea of CS divergence can be generalized to $L_p$ space by taking use of the Holder inequality: $|\int p(x)q(x)dx|\leq (\int|p(x)|^a d(x))^{1/a} (\int |q(x)|^b d(x))^{1/b}$, in which $a,b>1$ and $\frac{1}{a}+\frac{1}{b}=1$.

Hence, a generalized way to quantify the closeness between $p$ and $q$ can be expressed as: $D = -\log\left( \frac{|\int p(x)q(x)dx|}{(\int|p(x)|^a d(x))^{1/a} (\int |q(x)|^b d(x))^{1/b}} \right)$.

CS divergence is obtained by setting $a=b=2$. The generalized divergence can also be estimated by kernel density estimator.

[Q2 generalization of Theorem 1]

Thanks for the good suggestion. We generalized our Theorem 1 to any two square-integral densities $p$ and $q$. Please refer to Appendix B.5 and Proposition 5 in the revised submission.

[Q3 the same constant order of $I_{CS}(x;t)$ and $I(x;t)$]

Yes, the same constant order. By definition, $I_{CS}(x;t) = D_{CS}(p(x;t);p(x)p(t))$ and $I(x;t) = D_{KL}(p(x;t);p(x)p(t))$, our justifications in Appendix B.5 suggest that for any two square-integral densities $p$ and $q$, $D_{\text{CS}}(p;q) \lesssim D_{\text{KL}}(p;q)$. There is no other scaling factors in this inequality. By letting $p=p(x;t)$ and $q=p(x)p(t)$, we obtain $I_{CS}(x;t) \lesssim I(x;t)$.

[On the misalignment between the claims of contribution]

[P1: How accurate is the claim?]

Thanks for raising such a good question. More precisely, we feel that the independence to variational approximation and distributional assumption can be attributed to the estimation of CS divergence in a non-parametric way. Per your request, we analyzed how to non-parametrically implement KL divergence-based IB approach, from which we found the model is hard to train (see our reply to your [Q1]). We also point out that one possible reason why the non-parametric estimation of CS divergence successes is maybe because CS divergence relax the constraint on the support of distributions, in which we seldom has the situation $p(x) \log(\frac{p(x)}{0}) = \infty$.

[P2: Extending Theorem 1 to general distributions without Gaussian assumptions]

Theorem 1 can be generalized to arbitrary square-integral distributions without Gaussian assumptions. Please refer to Proposition 5 in the Appendix B.5 of the revised submission.

[P3: Strengthen the justification on adversarial robustness]

Thanks for your suggestion. We added Proposition 4 in the revised submission, in which we bound demonstrate that the probability that CS-QM, i.e., $I_{CS}(x;t)$ bounds the expected value of $|h_\theta(x+\delta)-h_\theta(x)|$ is at least $0.5$ and can reach to nearly $1$ (given sufficient samples). We also provide the assumptions under which this proposition holds.

[A lack of theoretical guarantees on the properties]

[P1&P2 The bias from non-parametric estimator and variational approximation]

First, thanks for raising such a good question. The bias and asymptotic property of our estimator has been discussed in Appendix B.4. The bias of existing variational approximation is also discussed in Appendix B.2 in the revised submission. In general, our argument is that previous methods either have much more obvious bias than ours or similar bias to us. We also provide a toy example in Appendix C.5 to show CS divergence estimated non-parametrically performs better than the variational KL divergence.

[P3 Justification of correlation with generalization performance in real-world data]

Thanks for the good suggestion. We justified the positive correlation on California housing data and also trained nearly 100 models with various depth, width, or mini-batch size. The results are demonstrated in Table 1 and Fig. 12. The experimental details are described in Appendix C.1.

[Concerns on regarding the reliability of results]

We added the standard deviation in Tables 8 and 9 in Appendix C.4. We explained in reply to your [Q5], there is no bug in the implementation. We also explained how hyperparameters in Table 4 are obtained in reply to your [Q7]. Please also refer to Appendix C.2.

Overall, many thanks for your careful reading of our submission. We hope that our reply and the revised submission can address your concerns.

[Q3: Fundamental difference with respect to KDE estimator in [3]].

Thanks for pointing out this good question; also sorry for missing this relevant discussion in the submission. Let us take the (marginal) CS divergence as an example, which can be expressed as: $D_{cs}(p;q)=\log\int p(x)^2dx+\log\int q(x)^2dx-2\log\int p(x)q(x)dx$

We focus our discussion on the first term inside the log, i.e. $\int p^2(x)dx$.

In fact, we estimate $\int p^2(x)dx$ as $\mathbb{E}_p(p(x))$. which is then estimated by $\frac{1}{N} \sum_i^N \hat{p}(x_i)$. This estimator is also called the resubstitution estimator; whereas [3] uses the plug-in estimator that inserting KDE of $p(x)$ into $p^2(x)$. That is, they estimate $\int p^2(x)dx$ as $\int\hat{p}^2(x_i)$, and use the property that “the integral of the product of two Gaussians is exactly evaluated as the value of the Gaussian computed at the difference of the arguments and whose variance is the sum of the variances of the two original Gaussian functions”. To our knowledge, this is a unique property that only holds for Gaussian kernel.

Please also refer to our discussion in Appendix B.4, especially the sentences marked with brown.

Based on the above analysis, we summarize the following two big differences:

1. although the final expression of CS divergence is similar. Our resubstitution estimator can be proven consistent (see also refer to Appendix B.4 in the revised submission), such property is missing in [3].

2. Moreover, our method is more general that is not restricted to Gaussian kernel (we do not to use any property associated with kernel function, such as convolution of two Gaussian).

[Q4: asymptotic properties of our estimator]

Yes, the asymptotic property of our estimator can be guaranteed. Pease refer to our Appendix B.4.

We start our analysis on the basic CS divergence: $D_{cs}(p;q)=\log\int p(x)^2dx+\log\int q(x)^2dx-2\log\int p(x)q(x)dx$

Proof sketch: we first provide the asymptotic property of each term inside the “log” operator. We obtained that for $\int p^2(x)dx$, $\int q^2(x)dx$ and $\int p(x)q(x)dx$, our empirical estimator combined with KDE has the asymptotic mean integrated square error (AMISE) tends to zero when the kernel size $\sigma$ goes to zero and the number of samples $N$ goes to infinity with $N^2\sigma\rightarrow\infty$. That is, each term has a consistent estimator. We then show that, the “log” operator does not influence the asymptotic property.

[Q5: Why performances of existing methods match perfectly in Table 2?]

First, there is no bug in our code. The performances of all existing methods are the same when the compression ratio $r=0$, i.e., no compression term exists in the IB objective, i.e., $\beta=0$. In this case, all existing methods just minimize a MSE loss (without any regularization on $I(x;t)$) with the same network architecture. Hence, their performances are the same.

Our CS-IB outperforms consistently when $r=0$. This results also demonstrate the advantage of using $D_{CS} (p(y|x);q_\theta(\hat{y}|x))$ alone over the traditional MSE loss.

[Q6: the standard deviation of results]

We provide the full results which include standard deviation (over 5 independent runs) in Tables 8 and 9 in Appendix C.4. Due to page limit, we only demonstrate the mean value in main paper.

[Q7: How hyperparameters in Table 4 were selected?]

For each competing method, the hyperparameters are selected as follows. We first choose the default values of $\beta$ (the balance parameter to adjust the importance of the compression term) and the learning rate $lr$ mentioned in its original paper. Then, we select hyperparameters within a certain range of these default values that can achieve the best performance in terms of RMSE. For VIB, we search within the range of $\beta \in [10^{-3}, 10^{-5}]$. For NIB, square-NIB, and exp-NIB, we search within the range of $\beta \in [10^{-2}, 10^{-5}]$. For HSIC-bottleneck, we search $\beta \in [10^{-2}, 10^{-5}]$. For our CS-IB, we found that the best performance was always achieved with $\beta$ between $10^{-2}$ to $10^{-3}$. Based on this, we sweep the $\beta$ within this range and select the best one. The learning rate range for all methods was set as $[10^{-3}, 10^{-4}]$. We train all methods for $100$ epochs on all datasets except for $PM2.5$, which requires around $200$ epochs of training until converge. All the hyperparameter tuning experiments are conducted on the validation set.

First, we would like to appreciate your careful review about our submission, and raise some insightful questions and suggestions. In the rebuttal, we systematically analyze the asymptotic property of our CS divergence estimator (please refer to Appendix B.4). We also analyze the sources of bias of previous variational KL divergence-based deep information bottleneck approaches (please refer to Appendix B.2). The advantage and distinctions of CS divergence over traditional distance measures like KL divergence is further introduced in Appendix B.1.

Moreover, per your request, we strengthened our theoretical results. First, we relax Gaussian assmption on the Theorem 1, we show in the Appendix B.5 that the relationship between CS divergence and KL divergence still holds for two arbitrary square-integral densities. An empirical justification is also provided. Second, we add Proposition 4, showing that our CS-QMI, i.e., $I_{CS}(x;t)$ bounds the expected value of $|h_\theta(x+\delta)-h_\theta(x)|$ with a probability of nearly $1$ (given sufficient samples). We also provide the assumptions under which this proposition holds.

Below, please let us reply to your questions one by one. We will then address your other concerns also one by one.

[Q1: non-parametric estimator of KL divergence]

Yes, the reviewer is correct that the conventional information bottleneck (IB) objective based on the KL divergence (i.e., Eq.~(16)) can also be evaluated non-parametrically with kernel density estimator (KDE). $D_{KL} (p(y|x);q(\hat{y}|x)) + \beta I(x;t)$.

The above objective actually equals to (by the chain rule of KL divergence): $D_{KL} (p(y,x);q(\hat{y},x)) + \beta I(x;t)$

The first term can be estimated by KDE with: $D_{KL} (p(y,x);q(\hat{y},x)) = \frac{1}{N} \sum_{i=1}^N \log \left( \frac{p(x_i,y_i)}{q(x_i,\hat{y}_i)} \right)$

In a Gaussian noise channel, $t=f(x)+\epsilon$ and $\epsilon \sim N(0,\sigma^2I)$. Hence, the second term equals to: $I(x;t)=H(t)-H(t|x)=H(t)-H(\epsilon)=H(t)-d(\log \sqrt{2\pi}\sigma + 1/2)$, in which $H(t)$ can also be measured with KDE.

We discussed this baseline in our revised submission. Please refer to our discussions on Appendix B.3 (especially the last paragraphs marked with blue). Unfortunately, we found that KL divergence estimated in a non-parametric way is hard to convergence or achieve a compelling performance.

One possible reason is perhaps the support constraint on the KL divergence to obtain a finite value: there is no guarantee that support of $q_\theta$ is a subset of the support of $p$. Our CS divergence does not have this issue, it relaxes the constraint on the support (unless two supports have no overlap). Please also refer to our discussion on Appendix B.1, in which we also plot Fig. 5 to clarify this point.

[Q2: CS divergence is not a special case of the well-known Renyi’ divergence]

First, we would like to make a clarification that the CS divergence is not a special case of the well-known Renyi’s divergence developed by Renyi with the following expression: $D_\alpha(p;q) = \frac{1}{\alpha-1}\log\int p^\alpha q^{1-\alpha}du=\frac{1}{\alpha-1}\log\int p\left(\frac{q}{p}\right)^{1-\alpha}du=\frac{1}{\alpha-1}\log\mathbb{E}_p\left(\frac{q}{p}\right)^{1-\alpha}$

When $\alpha=2$, we obtain a new divergence $-\log(1-\chi^2(p;q))$ that is proportional to the $\chi^2$ divergence, rather than our CS divergence mentioned in the paper.

To our knowledge, there are two big ideas to define divergence: the integral probability metrics (IRM) aim to measure $p-q$, whereas the Renyi’s divergences with order $\alpha$ aim to model the ratio of $\frac{p}{q}$. The general idea of CS divergence is totally different: it neither horizontally computes $p-q$ nor vertically evaluates $\frac{p}{q}$. Rather, it quantifies the “distance” of two distributions by quantify the tightness (or gap) of an inequality associated with integral of densities (i.e., $\int pq$). That is, CS divergence does not belong to either IRM or the traditional Renyi’s divergence family. Please refer to our discussion on Appendix B.1 for more details.

Our CS divergence belongs to a special case of the generalized divergence defined by Lutwak [Lutwak et al., 2005], which uses a modification of the Holder inequality. Lutwak et al. still use the name of “Renyi” maybe because their divergence also provides a one-parameter generalization of the KL divergence. But Lutwak’s definition is totally different from the traditional Renyi’ divergence.

To our knowledge, the use of the above Holder divergence or CS divergence is very scarce in the community. However, this family of divergence indeed has many desirable properties (please also refer to our discussion in Appendix D). We choose the most popular one in this new family of divergence in our paper. We also revise some relevant descriptions to avoid confusion.

First, thanks for your positive comments. Please let us reply to your concerns one by one.

[Q1: minor performance improvement]

Yes, on some relative simpler data, our prediction performance improvement seems a bit marginal. However, we would like to emphasize the obvious performance improvement on the adversarial robustness in Table 3.

Moreover, in the rebuttal period, we additionally added a new application: to predict the age of patients based on their brain functional MRI (fMRI) data with a graph neural network (GNN), which, to the best of our knowledge, is seldom investigated (prevalent literatures always rely on structural MRI (sMRI) data with convolutional neural networks) and is a much more challenging task. Our new results in Appendix D also demonstrate an obvious performance gain. Moreover, the success of our method in GNN also showcase the broader applicability and generalizability of our estimator.

[Q2: The gap between CS divergence and KL divergence]

$D_{KL}(p;q)$ is easy to reach infinity if the support of $p$ is not a subset of the support of $q$, which means that $D_{KL}(p;q)$ only has finite values if $\text{supp}(p) \subseteq \text{supp}(q)$. In fact, even if the support of $p$ and $q$ has large overlap, but $\text{supp}(p) \not\subset \text{supp}(q)$ (this situation is common in practice), KL is still infinity. Given two square integral densities $p$ and $q$, its value is infinite for CS divergence if and only if there is no overlap on the supports, i.e., $\text{supp}(p) \cap \text{supp}(q) = \emptyset$. In this sense, the CS divergence actually relaxes the support constraint of the KL divergence. Please also refer to our Appendix B.1 for more details.

[Q3: Why Cauchy-Schwarz divergence]

Thanks for the good question.

On the one hand, in Appendix B.2, we show that the approximation error of variational approximation is hard or intractable to control, or obvious larger than ours. By contrast, our estimation error and other properties can be fully analyzed as shown in Appendix B.4. That is, our estimator provides better theoretical guarantee.

On the other hand, some intrinsic limitations of KL divergence, such as support constraint, also motivates us to use other divergence. However, we would also like to emphasize here that the CS divergence also has some other desirable mathematical properties over KL divergence (as has been discussed in Appendix D), which we feel should be interesting and helpful to community.

[Q4: Toy examples where CS divergence performs slightly better]

Thanks for the good question. We consider two toy examples in the revised submission in Appendix C.5. The first example demonstrates when the true underlying distribution is a mixture of Gaussian, our CS divergence optimization performs better than the reverse KL divergence with variational approximation. The second example demonstrates that CS-QMI (the Cauchy-Schwarz quadratic mutual information) has similar or slightly better ability to identify complex dependence patterns compared with popular approaches like HSIC and distance covariance (dCov).

First, thanks for your positive comments. Please let us reply to your concerns one by one.

[Q1 The computational complexity of traditional KL-divergence based IB method]

Thanks for the good question. For HSIC-bottleneck, the dependence between $x$ and $t$ is measured with HSIC, which also takes a computational complexity of $O(N^2)$. For nonlinear information bottleneck (NIB) and square-NIB and exp-NIB, $I(x;t)$ is upper bounded as: $I(x;t) \leq -\frac{1}{N} \sum_{i=1}^N \log \frac{1}{N} \sum_{j=1}^N \exp\left( -\frac{\|h_i-h_j\|_2^2}{2\sigma^2} \right)$. in with $h$ is the hidden activity corresponding to each input sample, i.e., $t$ equals $h$ plus Gaussian noise (to avoid infinite mutual information value). Hence, the complexity is also $O(N^2)$.

[Q2 The effect of $\sigma$ to IB performance]

We demonstrate in Fig. 19 in the Appendix C.4, showing that our estimators provide consistent performance gain in a reasonable range of kernel size.

Additionally, we also demonstrate in Appendix B.4 that the effects of $\sigma$ to the bias and variance of our CS divergence estimator. The conclusion is that estimation bias is shrinking at a rate of $O(\sigma^2)$, and the variance shrinks at rate $O(1/(N^2\sigma))$. So there is a trade-off between bias and variance, which makes sense.

[Q3 The tuning of hyperparameter $\sigma$]

We partially agree with this point. For the basic variational information bottleneck (VIB) based on KL divergence, it is indeed very friendly for optimization. However, we would like to emphasize here that, for advanced KL-based IB approaches such as NIB, square-NIB and HSIC-bottleneck, they all require a careful treatment on the choice of some hyperparameter similar to kernel size $\sigma$. Different to these approaches, we systematically analyzed how $\sigma$ affects the bias and variance of our estimator in Appendix B.4. To our knowledge, such analysis is missing in previous literature.

[Q4 Potential disadvantage of CS-IB]

Thanks for the good question. We would like to mention a limitation of CS divergence with respect to the KL divergence. That is, we did not find a dual representation for the CS divergence. The dual representation of KL divergence, such as the well-known Donsker-Varadhan representation, enables the use of neural networks to estimate the divergence itself, which significantly improves estimation accuracy in high-dimensional space. A notable example is the mutual information neural estimator (MINE). We leave it as a future work and mention this point in Appendix D.