Asymptotics of Alpha-Divergence Variational Inference Algorithms with Exponential Families

Thank you for your thorough review and honest feedback. Below, we provide detailed responses to your second question. The two weaknesses you pointed out are adressed in our global rebuttal. We hope that our responses will clarify the issues you raised.

What is the technical difficulty to derive a convergence result similar to Theorem 1 for the proposed algorithm?

It does not present any major difficulty and can be done by following the exact same steps as in the proof of Theorem 1. We would get a similar result, with a slightly slower rate however: the $\gamma$ in the convergence speed would be replaced by $\gamma V(\eta_*)$, and we have $V(\eta_*) \leq 1$ by Jensen’s inequality. It is easy to gain some intuition for why this is the case: under the mean parameterization, the update is changed from $(1 - \gamma) \mu + \gamma \mathcal R$ to $(1 - \gamma V) \mu + \gamma V \mathcal R$, so we are still performing a convex combination of $\mu$ and $\mathcal R$ but with a less aggressive coefficient, thus yielding a slower rate of convergence. This also provides a heuristic as to why UNB moves slower at the start of the experiments: when the variational parameter is poorly chosen, $V$ is closer to $0$, which caps the progress that can be made at each step when the only values allowed for $(\gamma_t)$ are between $0$ and $1$.

We are grateful for your careful review, we hope our responses have adequately addressed your questions and provided further insights.

Thank you for your meticulous review and kind feedback. Below, we provide detailed responses to each of your questions, aiming to further clarify our work and address any remaining uncertainties.

Diverse remarks on notation + L196 : Notation $\circ$ is undefined (is this a composition?)

Yes, $\circ$ denotes a composition. We thank you for carefully reading our paper and will make sure to correct these mistakes.

L69 : Shouldn't the posterior be proportional to the joint $p(x, y)$ rather than the marginal $p(y)$?

Indeed, we can take $p(y) = p_{\theta}(x, y)$. The dependency in $x$ is dropped for notational ease.

L80 : "Assuming the argmax is uniquely defined at each iteration" it is unclear whether this is a reasonable assumption.

As you point out, this assumption can be challenging in the general case, but when the variational family is an exponential family as in (H1), the argmax in (1) is uniquely defined by $(1 - \gamma) \partial A(\eta) + \gamma \mathcal R(\eta)$ if and only if this quantity is in $F$. This is why the choice of $\gamma$ is so critical: if $\mathcal R(\eta)$ is not in $F$, we must take $\gamma$ small enough to ensure that the next iterate will be valid. Under (H1), the set $F$ is open and convex, which implies two things. First, by convexity, if $\gamma_0 \in (0, 1]$ is a valid choice, then all $\gamma \in (0, \gamma_0]$ are also valid choices. Second, by openness, it is always possible to find such a $\gamma_0$.

L108 : Is it obvious that $\mathcal R = U / V$?

Note that $V$ is the normalizing constant that makes $q_{\eta}^{\alpha} p^{1-\alpha} / V$ a probability density function. By linearity of the integral, $U / V$ is the expectation of the measurable function $S$ under a distribution with density $q_{\eta}^{\alpha} p^{1-\alpha} / V$, i.e. $\check\varphi_{\eta}^{\alpha}$. This explains why $\mathcal R = U / V$.

[Sec. 5.] falls a little flat as it concludes "we both obtain a biased and unbiased algorithm" without making those algorithms explicit or obvious. Also, it is unclear how strong an assumption it is that the exponential family does not depend on $x$ as in $q = q(\cdot | x)$.

The assumption $q_{\eta}(\cdot | x) = q_{\eta}$ in Section 5 was used only to illustrate a specific point: when $p_{\theta}$ is fixed and the decoder does not depend on $x$ (i.e., $p_{\theta} = p$ and we only update $\eta$), the exact gradient of the VR bound equals (up to a negative factor) $\mathcal R(\eta) - \partial A(\eta)$. This was meant to show that, in this particular scenario, an iteration of the algorithm discussed in Section 3 corresponds to a gradient ascent step on the VR bound. This assumption is not representative of the VAEs' training process, which uses the gradient estimators (17) and (18) with an Adam optimizer. Specifically, estimator (18) is used for VR and estimator (17) for UB. We will reorganize this section to clarify these points.

Thank you once again for your detailed review and supportive comments. We hope that our responses have satisfactorily addressed your questions and concerns.

Thank you again for your positive feedback. Regarding the notation, we understand your concern about potential confusion. We will emphasize the fact that the data $x$ is fixed throughout the optimization process, and will rename the prior $p_0$ to avoid any ambiguity and improve the accessibility of the paper.

We are grateful for your careful reading and constructive comments on our submission. Below, we respond to each of your questions, hoping to clarify any uncertainties.

The paper lacks clarity and is difficult to read.

We are sorry to hear that despite our best efforts to make the paper and the notation as clear and readable as possible, you encountered some difficulty in reading our work. We acknowledge that the notation in the paper is dense and will do our best to improve the overall readability by following your advice.

In Section 5, I was puzzled by the assumption that $q_{\eta}(\cdot | x) = q_{\eta}$, that is the factor does not depend on $x$. Does this mean that the encoder maps each image to the same latent representation?

The assumption $q_{\eta}(\cdot | x) = q_{\eta}$ in Section 5 was used only to illustrate a specific point: when $p_{\theta}$ is fixed and the decoder does not depend on $x$ (i.e., $p_{\theta} = p$ and we only update $\eta$), the exact gradient of the VR bound equals (up to a negative factor) $\mathcal R(\eta) - \partial A(\eta)$. This was meant to show that, in this particular scenario, an iteration of the algorithm discussed in Section 3 corresponds to a gradient ascent step on the VR bound. This assumption is not representative of the VAEs' training process, which uses the gradient estimators (17) and (18) with an Adam optimizer. Specifically, estimator (18) is used for VR and estimator (17) for UB. We will reorganize this section to clarify these points.

Figure 2: font size should be font size of text. Could the box plots be replaced with trajectories and shaded intervals?

Like you suggest, we originally intended for Figure 2 to show full trajectories with shaded intervals, but the result was cluttered and difficult to interpret due to the overlapping of five lines and their respective intervals. Additionally, the trajectories were quite noisy, although Polyak-Ruppert averaging could mitigate this specific issue. We will make further efforts to improve the clarity of the plot.

Define the IWAE and VAE baselines. Does this simply correspond to minimizing KL(q||p) and KL(p||q)?

Yes, IWAE is obtained when $\alpha = 0$ and VAE corresponds to $\alpha = 1$.

We extend our gratitude for your careful review and insightful feedback, and we hope our responses have effectively addressed your questions and concerns. We acknowledge the issues you noticed in the writing of our paper and will make sure to fix them.

We acknowledge that Section 5 lacks clarity and appreciate the opportunity to resolve this issue. In what follows, we expand on the reasoning behind this part of the paper.

Consider an algorithm whose updates are of the form $\mu_{t+1} = \mu_t +\gamma_t \mathrm M_th(\mu_t)$, where $\mathrm M_t$ is a positive definite matrix and $h:F\to F$. Usually, $h$ is the gradient of some function $H$ that we want to maximize.

In our case, it is equivalent to consider $\mathrm M_t = \mathrm{Id}$ and $h=\partial_{\eta} H$, or $\mathrm M_t=\mathrm F(\mu_t)$ and $h=\partial_{\mu} H$, where $\mathrm F(\mu)$ is the Fisher Information Matrix of the model at $\mu$,. This arises directly from the chain rule, along with the identities $\mu = \partial A(\eta)$ and $\partial_{\eta}\mu=\mathrm F(\mu)$ recalled in Appendix A.1. Indeed, they imply $\partial_{\eta} H(\mu)=\mathrm F(\mu)\partial_{\mu} H(\mu)$.

Since the gradients of interest are easier to compute w.r.t. $\eta$, we will choose the first option. Specifically, we have $\partial_{\eta}q_{\eta}(y)=\left(S(y) - \mu\right) q_{\eta}(y)$, which leads to $\partial_{\eta} V(\mu)=\alpha \left(U(\mu) - \mu V(\mu)\right)$.

For the unbiased algorithm $(5)$, we set $h(\mu)=\frac{U(\mu)}{V(\mu)} - \mu$, which we recognize as the gradient of $H:\mu\mapsto\frac{1}{\alpha}\log V(\mu)$.

Multiplying $H$ by the positive constant $\frac{\alpha}{1-\alpha}$ and writing it in integral form, we find that $\frac{\alpha}{1 - \alpha} H(\mu)=\frac{1}{1-\alpha}\log \mathbb E_{y\sim q_{\eta}}\left[\left(\frac{p(y)}{q_{\eta}(y)}\right)^{1-\alpha}\right]$ (remember that $\mu$ and $\eta$ are related by the one-to-one mapping $\partial A$).

In the more general case where we have $p_{\theta}(\cdot, x)$ instead of $p$ and $q_{\eta}(\cdot|x)$ instead of $q_{\eta}$, this is the Variational Rényi (VR) bound: $\mathcal L_{\alpha}^{\mathrm R}(p_{\theta}, q_{\eta}, x)=\frac{1}{1-\alpha}\log \mathbb E_{y\sim q_{\eta}(\cdot | x)}\left[\left(\frac{p_{\theta}(x, y)}{q_{\eta}(y| x)}\right)^{1-\alpha}\right]$.

In other words, the exact biased algorithm is a particular case of gradient ascent on the VR bound.

By doing a similar work for the unbiased algorithm, where $h(\mu_t)=U(\mu_t) - \mu_tV(\mu_t)$, we get $H(\mu)=\frac{1}{\alpha} V(\mu)=\frac{1}{\alpha} \mathbb E_{y\sim q_{\eta}}\left[\left(\frac{p(y)}{q_{\eta}(y)}\right)^{1-\alpha}\right]$. Hence, the unbiased algorithm is a particular case of gradient ascent on a new bound given by $\mathcal L_{\alpha}^{\mathrm G}(p_{\theta}, q_{\eta}, x)=\frac{1}{1-\alpha} \mathbb E_{y\sim q_{\eta}(\cdot | x)}\left[\left(\frac{p_{\theta}(x, y)}{q_{\eta}(y| x)}\right)^{1-\alpha}\right]$.

Using the REINFORCE gradient and backpropagation, we train VAEs to maximize these bounds. The respective methods are referred to as VR and UB in the paper.

We hope this reformulation of the proof will help clarify any ambiguities. Should you have any further questions or require additional explanations, please let us know. Thank you again for your time and constructive comments.

Thank you for engaging. However these additional details do not provide the clarification I'm looking for. Namely, why do we need to study the case where $q_\eta (\cdot \mid {\bf x}) = q_\eta$ in order to derive a learning algorithm?

In the above derivation, I do not see where in the calculation of $\frac{\alpha}{1 - \alpha} H(\mu)$ we exploit the fact that $p(\cdot, x) = p$ and $q_\eta = q(x)$. And once we derive $H$ in this particular case, why can we backtrack to the more general case?

I believe the writing should make it clear why this detour is necessary, how is the simpler case exploited, and why can we then go back to the more general case. Would it not be possible to directly work in the more general case?

Using the REINFORCE gradient and backpropagation, we train VAEs to maximize these bounds.

The authors have not defined "reinforce" gradient, could you define it?

Thank you for your feedback. We appreciate the opportunity to clarify our approach.

In the paper, we distinguish between two contexts of variational inference: the traditional framework discussed in sections 3 and 4, and the more complex VAE training setup.

1. In the traditional setting, we want to approximate the true posterior $p_{\theta}(\cdot|x)$ by learning the parameters of a single distribution $q_{\eta}(\cdot|x)$. Since the data $x$ remains fixed throughout the entire training process, we drop all dependencies in $x$: we suppose that $p_{\theta}(\cdot|x)$ is known up to a constant (through the function $p$), and simply we denote $q_{\eta}$. This does not mean that the learnt parameter isn't dependent on the data. Rather, we are optimizing a parametric distribution over the entire dataset, e.g., coefficients in Bayesian logistic regression.

2. In contrast, VAEs involve optimizing both an encoder $q_{\eta}(\cdot|x)$ and a decoder $p_{\theta}(x|\cdot)$ simultaneously. One of the main goals with VAEs is to learn a meaningful latent representation of the data, so each data point must be effectively mapped to the latent space, hence the particular conditioning on $x$. Unlike the traditional setup, where we learn a single distribution, VAEs require joint optimization of two networks. This complexity necessitates adjusting the approach.

The derivation in our paper shows that the algorithms MAX and UNB perform gradient ascent on specific variational bounds, namely the VR bound and $\mathcal L_{\alpha}^{\mathrm{G}}$.

The VR bound, for example, is defined as $\mathcal L_{\alpha}^{\mathrm R}(p_{\theta}, q_{\eta}, x)=\frac{1}{1-\alpha}\log \mathbb E_{y\sim q_{\eta}(\cdot|x)}\left[\left(\frac{p_{\theta}(x, y)}{q_{\eta}(y|x)}\right)^{1-\alpha}\right]$.

If we are in setup 1, the only thing we can act on is the parameter $\eta$. Differentiating the previous formula w.r.t. $\eta$ shows that MAX is a form of gradient ascent on $\mathcal L_{\alpha}^{\mathrm R}$.

For VAEs, we need to optimize both the encoder and decoder networks, and we are not changing $\eta$ and $\theta$ but the weights of the networks, of which the parameters are functions (say $\eta=\eta(W)$ and $\theta=\theta(W)$, where $W$ denotes the weights of both networks). To handle this new task, we use the reparameterization trick for differentiable optimization and the REINFORCE gradient, i.e., we write $\partial_{W}f(x, y, W)=f(x, y, W)\cdot \partial_{W}(\log f(x, y, W))$ with $f = \left(\frac{p_{\theta(W)}(x, y)}{q_{\eta(W)}(y|x)}\right)^{1-\alpha}$ and $\partial_{W}$ denotes the derivative with respect to the networks' weights. If further details are needed, we can provide a more comprehensive derivation.

What ties the algorithms from sections 3 and 4, and the VAE optimization process we propose together is the fact that they minimize the alpha-divergence by maximizing the same variational bounds. The procedures that lead to maximizing the bounds, however, are very different.

We hope this clarifies the fact that these two frameworks are different and thus need be treated adequately. Thank you for your consideration.

We appreciate your patience and apologize for misunderstanding your point of confusion, which is absolutely valid. This specific sentence in the paper is indeed wrong, due to a poor formulation of the underlying idea.

Throughout the entire paper, the function $q_{\eta}$ always depends on $x$ and can be written $q_{\eta}(\cdot|x)$. For simplicity and to reduce notation clutter, we omit this dependency until Section 5.

In Section 5, the function denoted $q_{\eta}(\cdot|x)$ is not the density of a distribution in an exponential family parameterized by $\eta$. Instead, it could be expressed $\tilde q_{f_{\eta}(x)}(\cdot)$, where $f_{\eta}$ is a neural network. When we wrote "the parameter $\theta$ is not updated, $q_{\eta}(\cdot|x)$ belongs to an exponential family and $q_{\eta}(\cdot|x)=q_{\eta}$ and does not depend on $x$", we meant that we were temporarily reverting to the "traditional VI" setting. The goal there was to explain why MAX performs gradient ascent on the VR bound when we are in the traditional setting and $q_{\eta}(\cdot|x)$ is the density of an exponential distribution.

The exact reasoning that led to the VAE training algorithms is the one described in the comments above. To summarize, we start by noticing that in the setting of Sections 3 and 4, the algorithms MAX and UNB are gradient ascent procedures, respectively on the VR bound and $\mathcal L_{\alpha}^{\mathrm G}$. Training VAEs to maximize these bounds leads to the methods referred to as VR and UB.

We hope this explanation will satisfactorily address your concerns about this part of the paper. We will delete the problematic sentence and reorganize section 5 so that it better reflects the train of thought behind the VAE training methods we use. Thank you for your thorough review and for helping us ensure the clarity of our paper.

Thank you for your interest in our work and for giving us the opportunity to improve it. Below is a draft of the revised Section 5 you asked for. We hope that this revision will clarify the points you raised.

In this section, we explain how to transpose the algorithms presented in Section 4.1 to the training of Variational Auto-Encoders (VAEs) [22]. We start by showing that the exact versions of the biased and unbiased algorithms correspond to gradient ascent procedures on two different variational bounds. Let us first address the case of the biased algorithm. Recall that it writes $\displaystyle{\mu_{t+1}=\mu_t+\gamma_t\left[\frac{U(\mu_t)}{V(\mu_t)} - \mu_t\right]}$. For $\eta\in E$ and $\mu=\partial A(\eta)$, we define the Variational Rényi (VR) bound [26] by $\mathcal L_{\alpha}^{\mathrm R}(\mu)=\frac{1}{1-\alpha}\log\mathbb E_{y\sim q_{\eta}}\left[\left(\frac{p(y)}{q_{\eta}(y)}\right)^{1-\alpha}\right]=\frac{1}{1-\alpha}\log V(\mu).$

Noticing that $\partial_{\eta}V(\mu)=\alpha \left(U(\mu) - \mu V(\mu)\right)$, we can thus express an iteration of the biased algorithm as $\mu_{t+1}=\mu_t+\frac{\gamma_t\cdot\alpha}{1-\alpha}\partial_{\eta}\mathcal L_{\alpha}^{\mathrm R}(\mu_t)=\mu_t+\frac{\gamma_t\cdot\alpha}{1-\alpha}\mathrm{F}(\mu_t)\partial_{\mu}\mathcal L_{\alpha}^{\mathrm R}(\mu_t).$ Under (H1), the Fisher Information Matrix $\mathrm{F}(\mu_t)$ is positive definite, hence the biased algorithm is a gradient ascent procedure on the VR bound.

The unbiased algorithm, which writes $\mu_{t+1}=\mu_t+\gamma_t\left[U(\mu_t) - \mu V(\mu_t)\right]$, similarly amounts to performing gradient ascent on the variational bound $\mathcal L_{\alpha}^{\mathrm G}$ defined by $\mathcal L_{\alpha}^{\mathrm G}(\mu)=\frac{1}{1-\alpha} \mathbb E_{y\sim q_{\eta}}\left[\left(\frac{p(y)}{q_{\eta}(y)}\right)^{1-\alpha}\right]=\frac{1}{1-\alpha} V(\mu).$

In the context of VAEs [22], we learn both a probabilistic encoder $y\mapsto \tilde q_{f(x;\eta)}(y)$ and a probabilistic decoder $x\mapsto\tilde p_{g(y;\theta)}(x)$, where $\tilde q_{\eta}$ and $\tilde p_{\theta}$ are densities from families parameterized respectively by $\eta$ and $\theta$, while $x\mapsto f(x;\eta)$ and $y\mapsto g(y;\theta)$ are neural networks. For simplicity and to align with the usual notation for VAEs, we will denote $\tilde q_{f(x;\eta)}(\cdot)=q_{\eta}(\cdot | x)$ and $\tilde p_{g(y;\theta)}(\cdot)=p_{\theta}(\cdot | y)$. We will also use the shorthand $\phi=(\eta, \theta)$.

Since the biased and unbiased algorithms studied in the previous sections minimize the alpha-divergence by maximizing the variational bounds $\mathcal L_{\alpha}^{\mathrm R}$ and $\mathcal L_{\alpha}^{\mathrm G}$, we propose to train VAEs to maximize those same bounds. In this new setting, they write $\mathcal L_{\alpha}^{\mathrm R}(\eta, \theta, x)=\frac{1}{1-\alpha}\log\mathbb E_{y\sim q_{\eta}(\cdot|x)}\left[\left(\frac{p_{\theta}(x, y)}{q_{\eta}(y | x)}\right)^{1-\alpha}\right],$ $\mathcal L_{\alpha}^{\mathrm G}(\eta, \theta, x)=\frac{1}{1-\alpha}\mathbb E_{y\sim q_{\eta}(\cdot|x)}\left[\left(\frac{p_{\theta}(x, y)}{q_{\eta}(y | x)}\right)^{1-\alpha}\right].$

To update both the encoder and decoder simultaneously, we differentiate them with respect to $\phi$ using the reparameterization trick [22]. If $z\sim r(\cdot)$ and there exists a mapping $v$ such that $v(z; \eta, x)$ has the same distribution as $y$ when $y\sim q_{\eta}(\cdot | x)$, then we have

$\partial_{\phi}\mathcal L_{\alpha}^{\mathrm R}(\eta, \theta, x)=\mathbb E_{z\sim r(\cdot)}\left[\overline w_{\alpha}(z, \eta, x)\partial_{\phi}\left(\log\frac{p_{\theta}(x, v(z; \eta, x))}{q_{\eta}(v(z; \eta, x) | x)}\right)\right],$

$\partial_{\phi}\mathcal L_{\alpha}^{\mathrm G}(\eta, \theta, x)=\mathbb E_{z\sim r(\cdot)}\left[w_{\alpha}(z, \eta, x)\partial_{\phi}\left(\log\frac{p_{\theta}(x, v(z; \eta, x))}{q_{\eta}(v(z; \eta, x) | x)}\right)\right].$

where $\displaystyle{w_{\alpha}(z, \eta, x)=\left(\frac{p_{\theta}(x, v(z; \eta, x))}{q_{\eta}(v(z; \eta, x)|x)}\right)^{1-\alpha}}$ and $\displaystyle{\overline w_{\alpha}(z, \eta, x)=\frac{w_{\alpha}(z, \eta, x)}{\int_{\mathsf Y} {w}_{\alpha}(z', \eta, x)r(z')\nu(\mathrm d z')}}$.

To train VAEs, we simply plug batch estimates of these gradients into an optimizer like Adam. Notably, $\partial_{\phi}\mathcal L_{\alpha}^{\mathrm G}$ can be estimated unbiasedly, while estimators of $\partial_{\phi}\mathcal L_{\alpha}^{\mathrm R}$ are subject to bias. We will study the practical implications of this fact in Section 6.

Thank you for your thorough and constructive review of our paper. We greatly appreciate the time and effort you have dedicated to providing feedback. We are committed to addressing the concerns raised and improving our work. Below, we respond to each of your points in detail.

The basis of the theoretical analysis is mostly covered by [1].

While our work indeed builds on the algorithm proposed in [1], which had already established the monotonicity property and explicit update form for the exponential family, our contributions lie in proving the algorithm’s convergence at an asymptotically geometric rate and towards a minimizer of the alpha-divergence. We further consider this algorithm in the empirical setting and propose an unbiased version, for which we provide two convergence theorems.

For Assumption (C2), the paper states that “$\alpha \in (0, ½)$ supposes specific behaviors on the relative tails”. It would be illuminating to verify this assumption in experiments.

We have found that convergence can still occur even if assumption (C2) is not met. For example, when $p$ is the density of a Cauchy distribution and the variational family is Gaussian, we still observe convergence (Appendix C). So far, we have not encountered a simple case where convergence entirely fails. This suggests that assumption (C2) could be relaxed or even replaced in future studies.

It would be nice to add experiments like Bayesian logistic regression.

Due to space constraints, we had to make choices regarding which experiments to present in the paper. We opted for the toy Gaussian examples as we found them to be more insightful than Bayesian logistic regression or other classical applications of variational inference algorithms.

The empirical results lack significance.

We believe that the main takeaway from the toy gaussian experiments is that the MAX approach is highly robust against aggressive hyperparameter tuning strategies, and generally converges faster than other methods with theoretical guarantees.

Although the NAT method performs well when it converges, it has shown occasional instability and is costlier per-iteration.

In VAE experiments, the UB approach is not significantly better than the VR approach in some cases.

It is true that the UB approach does not always show a significant advantage over the VR approach in our VAE experiments. However, we believe that this study offers useful insights. As you pointed out, the choice of $\alpha$ makes a difference in the empirical results, which is quite exciting as the idea behind the use of alpha-divergences in variational inference is to overcome some limitations of the traditionally used KL divergence which can be problematic on certain datasets. We observe that a proper tuning of $\alpha$ can greatly improve the model’s performance, though the optimal value seems to depend on both the gradient estimator and the dataset. A thorough analysis of the phenomena at play is beyond the scope of this paper, but is certainly a compelling question for future research.

The bias issue is not analyzed theoretically because “biased gradient estimators hinder any theoretical study”.

Analyzing the MAX approach theoretically presents two main difficulties. First, if it converges, it converges to a minimizer of an approximation of the alpha-divergence, rather than the true alpha-divergence itself. Second, this approximation is hard to characterize, since it amounts to finding and analyzing a function whose gradient is $\widehat U / \widehat V$.

Is it possible to compare the biased and unbiased algorithms on a relatively simple example like gaussian distributions or two-mixture of gaussians?

In the experiments section, both the biased and unbiased algorithms are compared in the case of a Gaussian variational family with Gaussian, mixture of Gaussian and Cauchy targets (see also Appendix C). The biased algorithm is referred to as MAX, while the unbiased algorithm is UNB.

If assumption (H3) is not satisfied, do we still have geometric convergence? If not, what would be the convergence rate?

We need assumption (H3) to obtain the geometric rate of convergence. Lemma 2 (found at the beginning of Appendix A.1) establishes that, as long as the sequence $(\eta_t)$ remains bounded, at least one fixed point of the mapping $\mathcal M_{\gamma}$ will be a limit point of this sequence. However, evaluating whether this limit point is the definitive limit and determining the convergence rate without (H3) is beyond the scope of this study.

UB achieved relatively worse results for $\alpha \in (½, 1)$, which is the interval where assumption C2 could be satisfied. Could the authors provide a heuristic understanding of this phenomenon?

We believe that the phenomenon at play is quite complex and depends mostly on the dataset, since opposite behaviors are observed between CIFAR10 and CelebA. Moreover, as we already discussed, assumption (C2) might be subject to relaxation, thus it may not fully explain the deeper reasons behind the observed convergence behavior.

Does [the monotonicity property] contribute to the empirical algorithm analysis? Is it beneficial to solving the instability problems for large step sizes?

Empirically, we can not guarantee the monotonicity property unless an accept-reject step is integrated into the procedure. The main idea behind the analysis of the sample-based algorithm is to use the mean parameterization of the exponential family to write it as a Robbins-Monro procedure, hence we do not leverage the monotonicity property seen in the exact setting. However, the fact that the exact versions of MAX and UNB enjoy this property may explain why these two approaches appear to be quite stable and follow direct paths to minimizers. This behavior is illustrated in Figure 1, and similarly observed in the mixture of Gaussian and Cauchy cases. We can include additional figures in Appendix C to further exemplify this matter.

Thank you for your thorough review and encouraging feedback on our submission. We provide detailed responses to each of your questions below, hoping to clarify any uncertainties you may have had reading our paper.

I am unsure why conditions (H4) and (C0') are considered realistic and sensible.

The rationale behind assumption (H4) is that as we approach the optimal solution, we can afford to be progressively less conservative in our choice of $\gamma$, since $\mathcal R(\eta_t)$ should always remain sufficiently close to the set $F$, or even belong to that set. More specifically, one can show that (H4) holds if the set $\mathsf K_{\mathcal L_{\alpha}(\eta_0)}$ is compact, meaning that the set of $\eta$'s that verify $\mathcal L_{\alpha}(\eta) \leq L_{\alpha}(\eta_0)$ is bounded. Consequently, each decrease in the alpha-divergence allows for an increase in the highest permissible value of $\gamma$. We are open to including a concise discussion on this matter in the Appendix.

Regarding (C0’), it is essentially a stricter version of (C0) and involves a design choice left to the practitioner's discretion. It's important to note that (C0) and (H4) are incompatible, which partly explains why we do not get a geometric convergence rate in empirical scenarios.

What influence does the number of samples $K$ have in the sample-based algorithm?

Increasing the number of samples $K$ reduces the bias of the estimator used in the first sample-based algorithm (MAX). Notably, the estimator $\widehat U / \widehat V$ becomes asymptotically unbiased, meaning that as $K \to +\infty$, it converges to $\mathcal R$. This is illustrated in the toy Gaussian experiment, where a small sample size ($K = 10$) causes the MAX algorithm to converge to suboptimal parameters with higher alpha-divergence compared to unbiased algorithms. To highlight the impact of sample size on bias, we can include additional plots in Appendix C. Aside from this, a smaller sample size accelerates per-iteration computation but results in noisier estimates of the intractable integrals. This creates a tradeoff, a detailed exploration of which is beyond this paper's scope.

Is it important to maintain the same number of samples as the number of iterations increases?

Theoretically, maintaining the same number of samples across all iterations is useful to satisfy condition (iv) in the proof of Theorem 3, as it guarantees the continuity of the function $\mu \mapsto \mathrm{Cov}(G(\mu))$. However, in practice, one might opt to use fewer samples in the early iterations and increase the number later on to improve the accuracy of the final outcome.

Thank you again for your attentive review and encouraging feedback on our submission, we hope our responses have helped resolve any ambiguities.