$t^3$-Variational Autoencoder: Learning Heavy-tailed Data with Student's t and Power Divergence

We are extremely grateful for your positive and through review of our work! We have updated the paper following your suggestions.

Addressing Weaknesses & Questions

(1) The authors claim that their idea is extendable to hierarchical models. Although there is a theoretical justification for this in the appendix, the paper would benefit from corresponding experimental studies.

We have newly implemented the two-layer hierarchical model ($t^3$HVAE) and compared against the Gaussian HVAE on double resolution images. The results are presented in Section 4.2, Figure 5 and Table C3. We demonstrate that the increased hierarchical depth allows $t^3$HVAE to learn more sophisticated images with substantially higher clarity and sharper detail compared to Gaussian HVAE, further justifying the generality and effectiveness of our theoretical framework.

(2) Minor comments: Thank you for going through our paper in detail! We have edited the mistakes and generally improved the presentation of the text. Also, please see our responses to the other reviewers for some other points that we addressed.

Thank you for your through review and suggestions, and raising a number of important issues! We have updated the manuscript to reflect our discussion, add experiments and improve presentation. Below are our responses to weaknesses which we hope can dispel any remaining concerns.

Addressing Weaknesses

• Does the $\gamma$-power divergence between the two manifolds still lead to the ELBO?

Thank you for pointing out the issue. We agree this is an important distinction that needs to be clarified, and the following discussion has been added to the manuscript in Section 3.2.

Due to the different expression for $\gamma$-power divergence and also the approximation step, the $\gamma$-loss does not function as a bound for log-likelihood, although joint minimization will still generally try to fit $p_\theta(x)$ to the true data distribution, maximizing likelihood while also keeping posteriors close. Its precise effect is a balance between reconstruction and regularization, which is explained throughout Section 3.3.

We emphasize that the central philosophy of our framework is viewing the ELBO not as a lower bound of likelihood (which naturally leads to modifying the KL regularizer term as most VAE variants do), but as a divergence between joint distributions (which naturally leads to modifying the entire divergence as in Eq.19). The $\gamma$-power or log-$\gamma$ divergences are reasonable choices due to their information geometric properties and viability of closed-form computations with respect to t-distributions.

But why does this work? Our experiments demonstrated that our framework is very effective even though we are not explicitly optimizing for likelihood. This can be explained in part by [1] which shows that $\gamma$-type divergence minimizers are M-estimators and proves asymptotic efficiency and robustness under heavy contamination, lending support to our approach as a solid alternative to MLE methods for statistical inference. Furthermore, for high-dimensional datasets $\gamma$ is quite small and the effect of changing divergence is simply less important compared to the effect of the coupled t-distributions (in particular the prior) whose degrees of freedom $\nu$ is small.

• Given that both t3VAE and the Student-t VAE integrate the Student's t-distribution into the VAE framework, a more granular comparison in the related work section would be enlightening.

We felt that the Student-t VAE design is less suitable to highlight in the introduction for a couple of reasons. The model uses a separate univariate t-distribution for each dimension (there is some confusion regarding dimension in the paper as well), and experiments only focus on overcoming the training instability issue for low-dimensional datasets. Moreover, performance improvement on image datasets in our experiments was not as significant compared to other models such as $\beta$-VAE and Tilted VAE. Similarly, DE-VAE also showed relatively poor reconstuction quality and VAE-st suffered from a range of inconsistencies and implementation issues. Hence we aimed for a broad overview of all the t-based VAEs rather than focusing on one single benchmark model.

• I noticed that the hierarchical variants of the model weren't evaluated in the experimental section. Such an evaluation might provide insights into the model's performance in different configurations.

We have newly implemented the two-layer hierarchical model ($t^3$HVAE) and compared against the Gaussian HVAE on double resolution images. The results are presented in Section 4.2, Figure 5 and Table C3. We demonstrate that the increased hierarchical depth allows $t^3$HVAE to learn more sophisticated images with substantially higher clarity and sharper detail compared to Gaussian HVAE, further justifying the generality and effectiveness of our theoretical framework.

• Incorporating the Student's t-distributions inherently increases the model's complexity. This can potentially introduce challenges in training and optimization.

Unlike models such as DE-VAE or VAE-st which require numerical integration to calculate the KL divergence between t-distributions for every data point, the explicit form of the $\gamma$-loss does not create any significant computational bottlenecks or new instability issues (e.g. division by zero variance, exploding gradients) compared to the original ELBO of Gaussian VAE. The same holds for the t-based reparametrization trick or multivariate t-distribution sampling process as detailed in Appendix C.1, requiring very minor additional computation. We also observed virtually equal runtimes in practice. Nonetheless we agree that this discussion is needed in the paper, and have merged into the final paragraph in Section 4 to give a general analysis on training cost and hyperparameter selection.

[1] H. Fujisawa, S. Eguchi, 2008. Robust parameter estimation with a small bias against heavy contamination. J. Multivar. Anal. 99(9).

Thank you for your positive review and very through advice which helped us greatly to improve our paper! We also report that we have additionally implemented the hierarchical $t^3$HVAE and added experimental results. Below are our responses to questions (merged with weaknesses).

• Is the gamma-power divergence not the same as the Renyi divergence (up to constants)? How do they differ?

In fact, the following four divergences have similar definitions.

$\alpha$-divergence $\frac{1}{\alpha(\alpha-1)}\left(\int p^\alpha q^{1-\alpha}-1\right)$

Renyi divergence $\frac{1}{\alpha-1}\log\int p^\alpha q^{1-\alpha}$

$\gamma$-power divergence $-\int p\left(\frac{q}{\lVert q\rVert_{1+\gamma}}\right)^\gamma +\lVert p\rVert_{1+\gamma}$

log-$\gamma$ divergence $\frac{1}{\gamma(\gamma+1)}\log\int p^{\gamma+1} - \frac{1}{\gamma}\log\int pq^\gamma+\frac{1}{\gamma+1}\log\int q^{\gamma+1}$

Despite the apparent similarities, they each possess rich geometric properties. The $\alpha$-geometry was originally studied by Amari and has analogous properties to KL divergence on the space of positive measures (being the intersection of the $f$- and Bregman families) but has a more complex dual connection when restricted to probability measures. The $\gamma$-geometry, studied by Eguchi, has a simple flatness structure where exponential families are replaced by power families and mixture geodesics are the same. Also, from a statistical perspective, $\alpha$ tunes between mass-covering or mode-seeking behavior, while $\gamma$ tunes between robustness and outlier-focused behavior [1].

In the initial stages of our research, we considered all the above as well as robust $\beta$ and other $f$-divergences as replacements for KL, but decided upon $\gamma$-power divergence due to flatness properties and viability of closed-form computations. The $\alpha$ and $\beta$ type divergences have also been studied for variational inference in previous works [2,3] (although not with our joint optimization framework), while $\gamma$-divergence methods are quite new.

• Can you provide an algorithm environment to show how this should be implemented at a glance?

Thank you for the helpful suggestion! We have added an implementation summary of our framework as Algorithm 1 in the main text.

• What do the MMD values themselves look like? These might be more useful to report than the p-values.

When the RBF kernel $k(x,x') = \exp\left(-\frac{(x-x')^2}{2\sigma^2}\right)$ is implemented for the MMD bootstrap test, the hyperparameter $\sigma^2$ is automatically chosen as the median of $\{(x_i - x_j')^2\}_{i,j}$. In practice, $\sigma^2$ fluctuates between around 3 to 5. In this case, directly comparing MMD values may not be meaningful. The MMD values are also somewhat randomly distributed around -0.01~0.01 and not particularly informative on their own. Hence a bootstrap test must be run to obtain meaningful statistics, which is why we only report $p$-values.

• Since you know the densities explicitly in the synthetic tests, you could use KSD instead of MMD, as this will have much improved statistical power. Have you tried this?

Due to increased power and reduction of reliable samples due to tail truncation, KSD tends to reject all models including $t^3$VAE. Of course, no model can exactly approximate the true distribution even with infinitely many samples, which will lead to rejection for any sufficiently powerful test. Instead our goal was to be able to clearly delineate relative tail generation performance of $t^3$VAE and other models, which is well achieved by MMD testing: we proved that $t^3$VAE is harder to distinguish from test data compared to other models. We also considered various other testing methods, but ultimately chose MMD due to applicability to 2 dimensions and simplicity & speed of implementation.

(continued below)