Inferring Principal Components in the Simplex with Multinomial Variational Autoencoders
Authors: James Morton
Justin Silverman
Gleb Tikhonov
Harri Lähdesmäki
Rich Bonneau
Venue:  (2021)
Keywords: Multinomial variational autoencoders
variational autoencoders
ILR transform
compositional PCA
probabilistic PCA
multinomial logistic normal

Abstract: Covariance estimation on high-dimensional data is a central challenge across multiple scientific disciplines. Sparse high-dimensional count data, frequently encountered in biological applications such as DNA sequencing and proteomics, are often well modeled using multinomial logistic normal models.  In many cases, these datasets are also compositional, presented item-wise as fractions of a normalized total, due to measurement and instrument constraints. In compositional settings, three key factors limit the ability of these models to estimate covariance: (1) the computational complexity of inverting high-dimensional covariance matrices, (2) the non-exchangeability introduced from the summation constraint on multinomial parameters, and (3) the irreducibility of the component multinomial logistic normal distribution that necessitates the use of parameter augmentation, or similar techniques, during inference. We show that a variational autoencoder augmented with a fast isometric log-ratio (ILR) transform can address these issues and accurately estimate principal components from multinomially logistic normal distributed data.
This model can be optimized on GPUs and modified to handle mini-batching, with the ability to scale across thousands of dimensions and thousands of samples.

Under review as a conference paper at ICLR 2021
INFERRING PRINCIPAL COMPONENTS IN THE SIMPLEX WITH MULTINOMIAL VARIATIONAL AUTOENCODERS
Anonymous authors Paper under double-blind review
ABSTRACT
Covariance estimation on high-dimensional data is a central challenge across multiple scientific disciplines. Sparse high-dimensional count data, frequently encountered in biological applications such as DNA sequencing and proteomics, are often well modeled using multinomial logistic normal models. In many cases, these datasets are also compositional, presented item- wise as fractions of normalized total, due to measurement and instrument constraints. In compositional settings, three key factors limit the ability of these models to estimate covariance: (1) the computational com- plexity of inverting high-dimensional covariance matrices, (2) the non-exchangeability introduced from the summation constraint on multinomial parameters, and (3) the irreducibility of the multi- nomial logistic normal distribution that necessitates the use of parameter augmentation, or similar techniques, during inference. Using real and synthetic data we show that variational autoencoder augmented with a fast isometric log-ratio (ILR) transform can address these issues and accurately es- timate principal components from multinomially logistic normal distributed data. This model can be optimized on GPUs and modified to handle mini-batching, with the ability to scale across thousands of dimensions and thousands of samples.
2 INTRODUCTION
Many scientific disciplines that collect survey data, such as economics, psychology, political science and the biological sciences routinely deal with compositional data, where only relative information can be measured. These datasets are often in the form of counts, where the total counts within sample are only indicative of the confidence of the measured proportions. The resulting proportions lie within a simplex and failing to account for the structure of this simplicial sample space can con- found the interpretation of the measurements. As result, there has been wide discussion across disparate disciplines 2 3 concerning the reproducibility crisis that has arisen from the misin- terpretation of compositional data. One of the obstacles the appropriate analysis of compositional data is the difficulty of efficiently estimating the latent parameters that lie in the simplex.
Accurately scaling probabilistic inference across high- nigh-dimensional count data is major outstand- ing challenge (5). This problem is apparent in the social sciences and is particularly pronounced in biological fields where datasets obtain observations tens of thousands of features across hun- dreds or millions of samples. One major computational bottleneck with Gaussian distributed data the inversion of d-dimensional covariance matrix that has a runtime of (dd) (6 7). As a re- sult, probabilistic covariance estimation for high-dimensional data is a computationally challenging problem.
Recent theoretical developments ( cementing the connection between Variational Autoencoders (VAEs) ( and Probabilistic Principal Components Analysis (PPCA) (10) holds much promise for enabling accurate, scalable, low-rank approximations of large covariance matrices. Variational au- toencoders were originally proposed as a generative model (9, but are now commonly deployed across scientific disciplines and have made contributions to single-cell RNA sequencing (11), mi- crobiome modeling (12), protein modeling processing (9) Following insights that connected regularized linear autoencoders and PCA (17) Lucas et al. 8 showed that carefully designed VAEs can recover the weights that are solved by PPCA. A computational advantage of VAES is that they do not require the inversion of a covariance matrix, and the resulting runtime (ndkIT) for n samples, d' dimensions, latent dimensions and T epochs. While it has been noted that VAEs may take tens of thousands of epochs estimate
(5), natural language processing (16) and image

---

Under review as a conference paper at ICLR 2021
the principal component (18), VAEs are easily parallelizable and can be accelerated with GPUs, pre- senting an attractive alternative to estimating principal components (17) and the resulting covariance matrix.
The connection between VAEs and PPCA is currently limited to Gaussian distributed data and not well-suited to a compositional setting. Showing that VAEs recover the correct principal com- ponents from count data is nontrivial due to the non-conjugacy issues between the logistic normal distribution and count distributions such as the multinomial distribution. Furthermore, the parame- ters of the multinomial distribution are compositional; they are constrained within the simplex and the resulting covariance matrix is singular and non-invertible 19). Aitchison (20) showed that PCA can be adapted to compositional data through the use of the center log-ratio (CLR) trans- form, which maintains isometry. However, this transformation not isomorphie, requiring that the resulting log-ratios sum to zero, and as result, CLR-transforme data will produce singular covariance matrix and rank-deficient principal components. has been shown that the isometric log-ratio (ILR) transform (2 satisfies both isomorphism and isometry and can handle this singu- larity issue (22 23) while enabling the estimation of full-rank principal components. Here, we show that VAEs augmented with the ILR transform can infer principal components learned from PPCA on multinomially distributed data, beginning to address these critical shortcomings.
RELATED WORK
In the microbiome literature, there have been number of methods (24 25 26 27 28 that have attempted to model ecological networks through the estimation of pairwise microbe correlations or pairwise inverse-covariance, where microbes are aggregated different taxonomical scales or 'taxa". Of these tools, only Flashweave can scale across more than thousands of taxa; however, it does this by avoiding the estimation of the covariance matrix. Methods that attempt estimate the covariance matrix can only handle on the order of a few thousand dimensions. Although there no widely accepted consensus definition of Multinomial PPCA in this context, being able to efficiently estimate the parameters of Multinomial PPCA would be highly useful for exploratory biological analysis. A number of studies have proposed using mixture modeling as a proxy for PCA |9|30) 30 31; however, these techniques depend either the Dirichlet distribution, whose covariance matrix is not flexible, or on stick bbrekknne which violates permutation invariance (32).
Lucas et al. ( has previously shown that the following two models can obtain the same maximum likelihood estimates of principal components W:
Probabilistic PCA
Linear VAE
p(z|z)= NW + n. I²) p(z) NO, I+)
p(z|z) N( u. o²,a (21xx) NVV(( ( p1) D)
Here, p(x[z) denotes the likelihood of observations æ Rd given the latent representation z RR p(z) denotes the prior on and q(2|z) denotes the estimated variational posterior distribution of 2 given an encoder parameterized by V and diagonal variances D. Both models estimate the same low dimensional representation of the data through 2, and learn the same factorization of the covariance matrix through W. While PPCA parameters are typically estimated through expectation maximization (1O), linear VAFs are optimized by maximizing the Evidence Lower Bound (ELBO) given by
logp(x) 2 E,(2(>>) [ogg(<//>]] - K [<({<=)))|pp((
For linear VAES with a Gaussian likelihood, the variational posterior distribution g( |x: can be shown to analytically agree with the posterior distribution p(=[æ) learned from PPCA (8). How- ever, deriving this connection for count-based likelihoods such as the multinomial distribution is complicated due to non-conjugacy issues (Appendix (A). This is major obstacle for many biologi- cal applications; multiple works have shown the merits of incorporating count distributions explicitly
into the model (33 34 35 36). Here, provide directions for overcoming this issue.
4 METHODS
First, we will redefine Multinomial PPCA with the ILR transform 21.. Then we will make the connection between Multinomial VAEs and Multinomial PPCA by leveraging insights from the
2

---

Under review as a conference paper at ICLR 2021
Collapse- Unnolllapee (CU) sampler (33). We will then derive algorithm to obtain the maximum posteriori (MAP) estimate for the VAE parameters.
4.1 PROBABILISTIC MULTINOMIAL PCA
PPCA can be extended to multinomially distributed data with the following generative model:
p(z|))) Mult( (Pn)
p(n]z)) = (W z, o Ia-1)
(2)
p(z) - IO,
(3)
Here W E Rd-1 represents the PCA loading matrix, o² is the variance, v E RRdxd- is a fixed contrast matrix whose columns sum to zero and & is the softmax transform. For a single sample, E Nd are the observed d -dimensional counts, E R= are the latent logits and z E R is the latent representation. The term @ (n)) distributed logistic normal, (m))) CN(W2, z, o² I). as shown by Aitchison (37). Furthermore, p(æ]2) yields multinomial logistic normal distribution, which is given by marginalizing out 1 the following expression:
MCN(x = a (p)()))))
This integral is not tractable; as result, this distribution does not have an analytically defined expec- tation, variance or probability density function. There have been multiple attempts to estimate the posterior distribution with MCMC (38) 33, (0)), but the complexity of this distribution requires a large number of samples, limiting the scalability of these methods. Variational methods have been developed to estimate the logistic normal distribution, but due conditional non-conjugacy. these methods often rely on approximations to the ELBO, further complicating estimation 411).
Recently, Silverman et al. 333 proposed use a Laplace approximation to estimate the parame- ters of the multinomial logistic normal posterior distribution. This approach relies on a two-stage optimization procedure on the factorized posterior distribution given by
p(7, =|x> x p(7)=))(z/())
(4)
If can be directly estimated, then conditional independence can be used to factorize the posterior distribution. Since the probability densities of the multinomial distribution and the normal distribu- tion are both log-convex functions, global optimum can be obtained for the multinomial logistic normal distribution (33) (Appendix B.3). Furthermore, the multinomial distribution does not in- troduce additional local optima for estimating Multinomial PCA. Given this, in addition to recent evidence that a PPCA MAP estimator can be obtained from regularized linear autoencoders (17), can design new algorithm to obtain a Multinomial PPCA MAP estimator.
4.2 THE ILR TRANSFORM ENFORCES IDENTIFIABLITY
The softmax function is a shift- iiftt-invariant function, which introduces an identifiability issue that has been addressed by the compositional data analysis community 42. In order remove the iden- tifiability issues, an isomorphism between the logits 77 and the multinomial parameters must be maintained. One commonly used solution is to use degenerate softmax, also known as the in- verse additive log-ratio (ALR) transform (I) (Appendix B.2).. Previous work has suggested that the isometric log-ratio transform (ILR) (21 (22 is more suitable for principal components analysis (Appendix B) The ILR and inverse ILR are given as follows:
ILR(x) " log x
ILR(x)) (())
where V Rdxd-1 is a basis such that ΦT Ia- and D = -,- flaxd A naive implementation of the ILR transform can be memory-intensive and computationally intensive for large d. However, any orthonormal basis can be used to parameterize the ILR transform and some of these bases can be represented by binary trees (43) 44 45). A binary tree can be used to represent T with Od d) elements, where the /th column vector of T is given as follows:
L V (0........... O, a a b O
(6)
k
a
(7)
3

---

Under review as a conference paper at ICLR 2021
where I indexes internal node in the tree with left children т, right children s. nodes to the left and nodes to the right t (46) (Figure S2). Due to rotation invariance, doesn't matter which tree is used to parameterize the ILR basis, but the choice of tree influence the runtime of the ILR transform.
a balanced binary tree is used, the memory requirements representing can be brought down from (dd)) to Oddlog d) and can reduce the matrix vector multplication runtime from ( to Old log d) (See Appendix B.I). This can speed the matrix-vector multiplication operations by an order of magnitude for datasets with more than ten thousand dimensions.
4.3 MULTINOMIAL VARIATIONAL AUTOENCODER ARCHITECTURE
Our full Multinomial VAE model is given follows:
p(a|n) (n) = Mult(o(WPn)
(8)
pl p/ dee) = NWW t u,
(9)
q(= 2z; Oene) (FF ( (log(o))))))))))))))))(()) (1)g D)
(10)
where Odec { o²) denotes the decoder parameters, Denc D} denotes the encoder parameters and μ E _d-1 is a bias parameter. Here, q==|2;; fene) denotes the variational posterior distribution of z given by the encoder represented as L-layer dense neural network with appro- priate activations. This encoder is directly used to evaluate p(n) 2; ddec). Furthermore, flat priors are assumed for all variables except z.
It is important to note potentially challenging modeling issues when designing the encoder. The ILR transform is not directly applicable count data, since log(O) is undefined. A common approach to this problem is to introduce a pseudocount before applying logarithm, which we will denote as log(æ)) (r) log( ( + 1). The choice of pseudocount is arbitrary and can introduce biases. To alleviate this issue, we introduce the deep encoder neural network highlighted in Equation I0, we expect that the universal approximation theorem would apply here (47 48) and that the accuracy of estimating the latent representation will improve with more complex neural networks. This is supported in our simulation benchmarks; more complex encoder architectures can better remove biases induced from the pscudocounts.
4.4 ALTERNATING LINI-BATCH OPTIMIZATION PROCEDURE
Given that our objective here is to obtain the MAP estimate of the VAE model parameters, the VAE parameters
{FL. W, D. o²) } can be obtained by estimating the global maximum of the posterior distribution. In the original CU sampler implementation, the parameter (71...., 7n is optimized across the entire dataset and then (7),- ....... 7m) fixed in order estimate the remaining parameters. For large studies, this can be memory demanding, since (71....,7m) 7m) E Baxd-1 alone can scale to millions of parameters.
To scale this estimation procedure to large high-dimensional datasets we have devised mini- batched alternating minimization procedure. For mini-batch X (a(1). ..... (())) of size b
Algorithm VAE Alternating Maximization Optimization
repeat for X (i) E X do
A < arg max
log ( H)
E bxbdd-1
"
0 arg max end for until convergence
log q( (
logp(Z)
o
and corresponding latent variables ( (2(2), ...., 2()) and (+)
1 (1)) the quan-
tities log ( H)) and log (0) are given by Equation 1 and The variational posterior
4

---

Under review as a conference paper at ICLR 2021
distribution of given by q H x ; 0) IT ((170)/2001) 0) be obtained by marginaliz-
ing out z in Equations 8 and Toas as follows:
q(n q(1 x
. (W 4 (log( u4) + M, WDW +
(11)
The prior logp(25 ) is evaluated, given log log () Mog.((Z' I&) and the latent encoding
means given by ( Fz ( ( (log(X (1) u)), where (6) is integrated out of Equation
RESULTS
5.1 MULTINOMIAL PPCA AND VAES AGREE IN SIMULATION BENCHMARKS
To determine if the proposed multinomial variational autoencoder can recover principal components, extended the benchmarks proposed in (8). Here, we benchmarked proposed analytical VAE
Analytical Mult. VAE
Stochastic Mult. VAE
Mult. PCA (MAP)
Mult. PCA (HMC)
a
b
C
Ground Truth Covariance
r=9.966
r=0.907
r=0.78
r=0.8
0.5
0.5
0.5
0.5
0.0
0.0
-0.5
-05
2
0.00
0.05
200
0.00
0.05
g Subspace distance
h
Correlation
Axis-alignment
Procrustes
0.100
0.100
0.100
0.100
Probability
0.050
0.050
0.050
0.050
0.025
0.025
0.025
0.025
0.8
0.9
0.6
0.7
0.10
0.11
0.2
0.4
0.6
VAE architecture
-wis-alignment
1.0
Subspace distance
stochastic-ILR
0.9
Procrustes
0.8
analytic-ILR
stochastic-ALR
0.8
0.4-
analytic-ALR
stochastic-log
0.2
analytic-log
50
100
50
100
Epochs
Epochs
Epochs
Dense simulation benchmarks. (a-d) Comparision of ground truth covariances (y-axis) and estimates from analytical VAE, stochastic VAE, MAP estimator and HMC performed in Stan (x-axes); all of these methods use the ILR transform and only the posterior means are shown in (d). (e-h) Correlation, alignnligmmnn, subspace distance and Procrustes metries between analytical VAEs, stochastic VAEs. MAP and HMC after convergence. (i-k) Comparison of log, ALR and ILR transforms with respect axis-alignment, subspace distance and Procrustes. The axis-alignment metric measures the angular differences between W and the ground truth principal components. The subspace distance is a measure agreement between the ground truth correlations and the correlations estimated from W. Procrustes measures the residual eπror after obtaining the best rotation and scaling factors to match W to the ground truth principal components. See Appendix C.5 for descriptions of the correlation, axis-alignment, dubspace-distance and Procrustes metrics.
Figure 1:
to the stochastic VAEs with a multinomial likelihood across multiple simulations. These method- ologies were compared against the ground truth covariances of the multinomial logistic normal, in addition to the MAP and Hamiltonian Monte Carlo (HMC) estimates obtained from Stan (49). The agreement between the ground truth principal components and the estimated principal components measured by axis-alignment, subspace distance, correlation and Procrustes (50) (see Appendix
C.4 for details).
5

---

Under review as a conference paper at ICLR 2021
When fitting against multinomial logistic normal distributed data with counts, all of the ILR- based methodologies can accurately estimate the covariance matrix up to scale (Figure [a-d). The analytical VAE MAP estimate and the posterior samples of HMC all have correlation close to 1, suggesting strong agreement between the ground truth covariance and the estimated covariances. the principal components
perfectly estimated, the axis-alignment, subspace distance and Pro- crustes metrics would all be close to zero. While the subspace distance and Procrustes are notably close to zero, the axis-alignment metric is above 0.5, which would suggest disagreement between the trouh and the estimated principal component axes. Given that the posterior distribution of the axis-alignment metries obtained from HMC overlaps with the analytical VAE, stochastic VAE and MAP estimates, there is evidence that reducing this metric to zero is inherently difficult for count data. However, both the subspace distance and Procrustes metrics approach zero, supporting our claim that principal components can be estimated data. Furthermore, our simulations (Figure Ti-k and Figure 21-h) suggest that among the log-ratio transforms benchmarked. only the ILR transform recover principal components, corroborating previous findings surrounding compositional PCA
scale and rotation on fully observed count
42 23 Dealing with sparse data presents addi-
Log-ratio Transform Benchmark
1.00
1.0
Axis-a aiis-alignment
distance
0.9
Correlation
0.8
Procrustes
0.90
0.85
0.80
0.6
0.5
500
1000
500
1000
500
1000
500 Epochs
1000
Epochs
Epochs
Epochs
Correlation
0.15
0.10
Axis allsamegnnt
0.996
0.998
Subspace distance
0.975
0.950
Procrustes
0.96
0.98
0.94
0,900
0.00
0.992
0.875
0.92
500
1000
500
1000
500
1000
500
1000
Epochs
Epochs
Epochs
Epochs
stochastic-ILR
stochastic-ALR
stochastic-log
VAE architecture
analytic-ILR
analytic-ALR
analytic-log
Encoder Depth Benchmark
Axis-alise alignment
1.00
Subspace distance
1.0
Correlation
Procrustes
0.90
1000
500
1000
1000
Epochs
Epochs
Epochs
Epochs
Subspace distance
Axiss dits-aligmmen
0.998
0.975
0.98
Correlation
Procrustes
0.10
0.996
0.994
0.94
0.900
0.992
0.875
0
500
1000
500
1000
500
1000
500
1000
Epochs
Epochs
Epochs
Epochs
VAE architecture
stochastic 1-layer
stochastic 3-layer
stochastic 5-layer
analytic 1-layer
analytic 3-layer
analytic 5-layer
Figure 2: Comparison of different log-ratio transforms (a-f) and encoder depths (g-l) in sparse sim- ulations. Two datasets with 200 and 5000 input dimensions are simulated for benchmarks (a-d, i-l) and (e-h, m-p) respectively. Log-ratio transform benchmarks investigate different log-ratio trans- forms using a single-layer encoder. The Encoder Depth benchmark investigates the impact of the number of encoder layers combined with the ILR transform.

---

Under review as a conference paper at ICLR 2021
tional challenges; the zeros in count data indicative of missing data, which can further complicate the estimation of principal components. Our hypothesis that boosting the complexity the encoder architecture would help alleviate issues with missing data is supported by benchmarks shown in Fig- ure -- In both of the sparse datasets, none of the methods were able to achieve optimal accuracy any of the benchmark metrics, but there is clear advantage of utilizing multilayer encoders compared to single-layer encoders.
Across the simulation benchmarks, the analytical and stochastic VAEs have comparable perfor- mance, with discrepancies highlighted by the correlation metric in the dense and sparse bench- marks. On the dense datasets, our proposed analytical VAE has better agreement with the ground truth covariance metric, whereas on the sparse datasets, it appears that the stochastic VAE has better agreement. It is difficult to explain the reason for the performance gap between our proposed ana- lytical VAE and the stochastic VAE due the analytical intractability of the Multinomial VAE ELBO. The challenge of accurately estimating an optimal H for each mini-batch could be one limiting factor affecting the performance of our proposed analytical VAE (Appendix B.4..
5.2 PRETRAINING MULTINOMIAL VAES ON VERY LARGE COMPENDIUM OF MOUSE MICROBIOME OBSERVATIONS
To showcase this on real microbiome datasets, 11k mouse fecal pellet samples with total of 8.8k features across 45 studies were analyzed. We trained 2 models to evaluate the differences between the stochastic VAE and the analytical VAE. We trained these models for 10k epochs; the models with the smallest validation error are reported here. The details behind the full training procedure given in Appendix C.3
Visualization of the learned decoder weights shown in Figure makes
a
Stochastic VAE
b
Analytical VAE
250
400000
Embedding dimension
200
Embedding dimension
30000
150
100
20000
50
0
10000
-50
-0
-100
Embedding dimension
Embedding dimension
d
0.07
diagonal off-diagonal
0.00020
0.06
0.05
0.00015
Density
Density
0.04
0.00010
0.03
0.02
0.00005
0.01
0.00
0.00000
-100
0
100
200
300
-10000
0 10000 20000 30000 40000
Figure 3: Comparison of the decoder inner product between the pretrained stochastic VAE and the pretrained analytical VAE, (a-b) Heatmap of WWW. (c-d) Distribution of diagonal elements and off-diagonal elements of WT W.
W.
it clear that the analytical VAE and, to lesser extent, the stochastic VAE are able learn orthog-

---

Under review as a conference paper at ICLR 2021
onal embeddings. since W W approximates diagonal matrix. The orthogonality of the decoder weights is more apparent in our proposed analytical VAE than in the stochastic VAE decoder; there is larger distinction between the diagonal elements and the off-diagonal elements. Unlike linear VAFs (8), the connection between the eigenvalues and the decoder weights is not as clear due to the asymmetry between the encoder and the decoder; scale identifiability of the decoder weights complicates the process of estimating the eigenvalues. This is apparent in Figure] heatmaps of the W inner products are on vastly different scales.
One of the advantages of employing a pretraining strategy with VAEs is that the pretrained models can enable one-shot learning (5I). By estimating a low-dimensional representation from a large unlabeled sigh-dimensional dataset, fewer labeled samples for training downstream classifiers We showcase this property in two classification benchmarks.
5.3 CLASSIFICATION AND CORRELATION BENCHMARKS
To construct benchmarks on real biological datasets, we analyzed 16S sequencing data obtained from study conducted by Tripathi et investigating hypoxia that used mouse fecal pellets (52)
and study conducted by Shalapour et al. investigating cancer effects mice (53)
Due to the connections between Multinomial VAEs and compositional PCA that, if our VAE is performing well, we will see that the row Euclidean distances of V E Raxk correspond to Lovell's un-normalized proportionality metric (54) 23) given by
20
22
we expect
(W); - (WW), x Var log Zi
(12)
2
IJ
where refers to the observed proportions and and j refer to the two features being compared. The proportionality metric has a simple interpretation; if this variance is small, that implies that the two features are highly co-occurring.
Here, only a relative relationship between the proportionality metric and the VAE decoder row dis- tances can be stated due to scale identifiability issues inherent in Multinomial PPCA. Proportionality has been shown to be a more reliable metric for compositional data compared to standard correla- tion metrics such as Pearson and Spearman single cell RNA sequencing and microbiome studies (,55 56 2), supporting the notion that our proportional comparison of our VAE decoder to Lovell's un-oormalized proportionality metric will provide sound perspective, as well means for inter- preting the weights of the VAE.
We compared the estimated VAE embeddings with the proportionality metric to determine if this relationship holds empirically. The pairwise decoder distances are compared to the pairwise propor- tionality metries of 300 selected features (Table and Figure 55). The learned VAE representations benchmarked using K-nearest neighbors (KNN) classification. KNN is applied to the learned VAE encodings across the two microbiome datasets to determine how well the classifiers can iden- tify the mice that were induced with hypoxia (52) and classify mice based on their experimental group (53). These classification models are compared against two baseline models, namely KNN trained on LDA topic proportions (31) and KNN trained on raw counts.
From the classification benchmarks illustrated in Table I we can see that there is a large margin between the VAE models and the baseline models across the measured metrics. Both VAEs had a higher Fl score and AUC, suggesting that the learned representations can better separate the ex- perimental groups. Part of this performance boost could be attributed to the differences between distance metrics. When the ILR transform is utilized, the Euclidean distance between the latent encodings approximates the Aitchison distance between proportions (57 which measures rel- ative differences instead of absolute differences. Our findings collaborate theoretical insights from compositional data analysis and empirical microbiome observations, where the Aitchison distance provided a substantial boost in KNN classification accuracy compared to other distance metrics (9). Furthermore, the average negative log-likelihood score (NLL) is lower for both of the VAE models compared to LDA, suggesting that the VAE models generalized better on held out data than LDA. This decrease in the predictive log- legllikllhoodd could be attributed to the increased flexibility of the covariance matrix in the logistic normal distribution compared to the Dirichlet distribution.
There are a couple of notable discrepancies between the two Multinomial VAE models. The stochas- tic VAE appears have superior classification performance and lower reconstruction emor com-
8

---

Under review as a conference paper at ICLR 2021
Dataset
Method
NLL
Fl score
AUC
Proportionality
Shalapour et al.
analytical- VAE
57750
0.700 = 0.021
0.910= 0.013
0.60
stochastic-VAn VAB
11615
0.732+ 0.023
0.891 0.015
0.30
LDA
169384
0.507 = 0.034
0.742 + 0.018
NA
raw
NA
0.406 = 0.039
0.686 + 0.018
NA
Tripathi al.
analytical- VAE
117789
0.867 = 0.005
0.873 + 0.004
0.62
stochastic-Va VAB
13815
0.926 0.004
0.928 0.003
0.28
LDA raw
132091
0.718 = 0.006
0.744 + 0.005
NA
NA
0.700 + 0.007
0.729 + 0.006
NA
Table 1: Classification and Correlation benchmarks. KNN classification was performed with k=5 for all representations. The average negative log- ligeliilldood (NLL) was evaluated samples that were held out during the training of the VAEs. 100 rounds of 10-fold cross validation were applied estimate average Fl score, AUC and standard errors. "Raw" indicates that transformation was applied to the raw counts before performing KNN classification. The agreement between the VAE decoder weights and Lovell's proportionality metric is measured according Pearson's r on the log- -egtrranformeed metrics.
pared to the analytical VAE. The exception to this is the HCC dataset (53), where the analytical VAE marginally outperforms the stochastic VAE in terms of AUC. However, the analytical VAE can learn more apparent orthogonal embeddings (Figure (3) and better agrees with Lovell's proportion- ality metric (Table Db. Only the analytical VAE was able recover the log-linear relations between Lovell's proportionality and the VAE embedding distances, suggesting that it can more accurately learn biologically relevant correlations (Figure S
CONCLUSION
Prior work aiming build frameworks for probabilistic estimation of covariance matrices from count data have been largely limited to conjugate priors due to their tractability. However, these choices can lead to models with lower explanatory power due the rigid structure of the resulting covariances. Further, the correct treatment of compositional data requires additional development in this context. For example, disregarding the simplicial sample space associated with composi- tional data and performing Pearson correlations on raw count data is common practice in scientific applications, but is a source of reproducibility issues (2). Due to the negative bias induced from the covariance the observed count data, the estimated covariances will not agree with the ground truth covariances in the system of interest, a fact noted by Pearson in 1897 (60).
Adapting the logistic normal distribution in place of conventional conjugate priors provides means to remedy the issue of inferring correlations on compositional data. The covariance matrix on ILR- transformed data can be interpreted using Lovell's proportionality metric, and can serve placement for pairwise correlations. Since log-ratios are scale- -cvle-invariant, the dependence on the total counts disappears, which is critical for assuring agreement between the relative quantities measured through count data and the absolute quantities in the system that aren't 't directly observable. While there is sound motivation has been limited due to the challenge of estimating these distributions. Here, we have provided a means to accurately and efficiently estimate these distributions and covariance matrices up to scale using the ILR transform. To this end, we have provided a proof- proof-of-coneeee -oneccc--ccptttt that Multinomial VAEs can learn the principal components obtained from Multinomial PPCA.
employ the logistic normal distribution in this context, its application
We have shown that fitting low-rank approximations using Multinomial VAEs can provide more explanatory representations than LDA while providing means to obtain interpretable correlations. These methods can be used in one-shot learning and transfer learning settings, requiring fewer la- beled samples for training and allowing for the use of pretrained models to extract features from smaller datasets. Given the vast number of scientific disciplines that collect compositional data, we anticipate that models such as these Multinomial VAEs will have a significant impact on many scientific applications.

---

Under review as a conference paper at ICLR 2021
CODE AVAILABILITY
All of our software and analyses can be found on Zenodo at http://doi.org/11.... 10..0//////00..0//// 1/ten/oo..2200000000
ACKNOWLEDGEMENTS
We want acknowledge Pytorch (61. Pytorch-Lightning (62) the Biom format (63. Matplodlib (64), Scipy (65), Numpy (66) and Stan 44) for providing the software foundation that this work was built upon.
REFERENCES
[1] John Aitchison. The statistical analysis of compositional data. Journal of the Royal Statistical Society: Series (Methodological), 1402):139-160, 1982.
[2] Gregory B Gloor, Jean M Macklaim, Vera Pawlowsky-Glan and Juan Egozcue. Microbiome datasets are compositional: and this is not optional. Frontiers in microbiology, 8:2224, 2017.
[3] Michael A Skinnider, Jordan W Squair, and Leonard J Foster. Evaluating measures of associ- ation for single-cell transcriptomics. Nature methods, 16(5):381-386, 2019.
[4] Donald A Jackson. Compositional data in community ecology: the paradigm or peril of pro- portions? Ecology, 78(3):929-944) 1997.
[5] David M Blei, Alp Kucukelbir, and Jon D McAuliffe. Variational inference: A review for statisticians. Journal of the American statistical Association, 112(518):859-8777 877, 2017.
[6] Steven C. Althoen and Renate Melaughlin. Gauss-jordan reduction: A brief history. The American mathematical monthly, 1402):130-142, 1987.
[7] Aravindh Krishnamoorthy and Deepak Menon. Matrix inversion using cholesky decomposi- tion. 2013 signal processing: Algorithms, architectures, arrangements, and applications (SPA), pages 70-72. IEEE, 2013.
[8] James Lucas, George Tucker, Roger Grosse, and Mohammad Norouzi. Don't Blame the ELBO! A Linear VAE Perspective Posterior Collapse. (NeurIPS), 2019.
[9] Diederik P Kingma and Max Welling. Auto-encoding variational bayes. arXiv preprint 131iv:1312.61 14, 2013.
[10] Michael E Tipping and Christopher M Bishop. Probabilistic principal component analysis. Journal of the Royal Statistical Society: Series (Statistical Methodology), 61(():611-62) 1999.
[11] Romain Lopez, Jeffrey Regier, Michael Cole, Michael I. Jordan, and Nir Yosef. Deep gen- erative modeling for single-cell transcriptomics. Nature Methods, 15(12):0033-1058 2018.
[12] Min Oh and Liqing Zhang. Deepmicro: deep representation learning for disease prediction based on microbiome data. Scientific reports, 10(1):1-9, 2020.
[13] Sam Sinai, Eric Kelsic, George M Church, and Martin A Nowak. Variational auto-encoding of protein sequences. arXiv preprint 17
2017.
[14] Xiaojie Guo, Sivani Tadepalli, Liang Zhao, and Amarda Shehu. Generating tertiary protein structures via an interpretative variational autoencoder. arXiv preprint arXhv:2004.07119, 2020.
[15] Alex Hawkins Hooker, Florence Depardieu, Sebastien Baur, Guillaume Couairon, Arthur Chen, and David Bikard. Generating functional protein variants with variational autoencoders. BioRxiv, 2020.
10

---

Under review as a conference paper at ICLR 2021
[16] Zichao Yang, Zhiting Hu, Ruslan Salakhutdinov, and Taylor Berg- Kirg-Kirpattrik. Improved variational autoencoders for text modeling using dilated convolutions. arXiv preprint arXiv:1702.08139, 2017.
[17] Daniel Kunin, Jonathan M Bloom, Aleksandrina Goeva, and Cotton Seed. Loss landscapes of regularized linear autoencoders, arXiv preprint arXiv:1901.08/ 2019.
[18] Xuchan Bao, James Lucas, Sushant Sachdeva, and Roger Grosse. Regularized linear autoen- coders recover the principal components, eventually, arXiv preprint arXiv:2007.0673/ 2020.
[19] Peter Filzmoser and Karel Hron. Correlation analysis for compositional data. Mathematical Geosciences, (1):909, 2009.
[20] John Aitchison. Principal component analysis of compositional data, Biometrika, 70(1):57-65, 1983.
[21] Juan José Egozcue, Vera Pawlowsky-Glabe Glòria Mateu-Figuern and Carles Barcelo-Vidde Isometric logratio transformations for compositional data analysis. Mathematical Geology, 35(3):279-300 2003.
[22] Peter Filzmoser, Karel Hron, and Clemens Reimann. Principal component analysis for com- positional data with outliers. Envirommetrics: The Official Journal of the International Envi- roumetrics Society, 20(6):621-632, 2009.
[23] Juan José Egozcue, Vera Saslowsky-Glahnee and Gregory B Gloor. Linear association com- positional data analysis. Austrian Journal of Statistics, 47(1):3-31, 2018.
[24] Zachary D Kurtz, Christian L Müller, Emily R Miraldi, Dan R Littman, Martin J Blaser, and Richard A Bonneau. Sparse and compositionally robust inference of microbial ecological networks. PLoS Comput Biol, 1165)ee0042266, 2015.
[25] Jonathan Friedman and Eric J Alm. Inferring correlation networks from genomic survey data. PLoS Comput Biol, 1021)77, 2012.
[26] Julien Chiquet, Mahendra Mariadassou, and Stéphane Robin. Variational inference for sparse network reconstruction from count data. pages 1-30, 2018.
[27] Huaying Fang, Chengcheng Huang, Hongyu Zhao, and Minghua Deng. geoda: conditional dependence network inference for compositional data. Journal of Computational Biology, 24(7):699-708, 2017.
[28] Janko Tackmann, João Frederico Matias Rodrigues, and Christian von Mering. Rapid infer- ence of direct interactions in large-scale ecological networks from heterogeneous microbial sequencing data. Cell systems, 93)::86--266, 2019.
[29] Wray L Buntine and Aleks Jakulin. Applying discrete pea in data analysis. arXiv preprint /207.4/ 2012.
[30] Wray Buntine. Variational extensions to em and multinomial pea. In European Conference on Machine Learning, pages 23-34. Springer, 2002.
[31] David M Blei, Andrew Y Ng. and Michael Jordan. Latent dirichlet allocation. Journal of machine Learning research, 3(Jan)993-1022, 2003.
[32] Scott W. Linderman, Matthew J. Johnson, and Ryan P. Adams. Dependent multinomial models made easy: Stick breaking with the Pólya-gamma augmentation. Advances in Neural Informa- tion Processing Systems, 2015 3-1aaaaaa:3446-3444 2015.
[33] Justin D. Silverman, Kimberly Roche, Zachary C. Holmes, Lawrence A. David, and Sayan Mukherjee. Bayesian Multinomial Logistic Normal Models through Marginally Latent Matrix- Processes. pages 1-39,2019 2019.
[34] Travis E Gibson and Georg K Gerber. Robust and scalable models of microbiome dynamics. arXiv preprint arXiv: 1805.945) 2018.

---

Under review as a conference paper at ICLR 2021
[35] Tarmo Aijö, Christian L Müller, and Richard Bonneau. Temporal probabilistic modeling of bacterial compositions derived from 16s mma sequencing. Bioinformatics, 34(3)377--880, 2018.
[36] Tyler A Joseph, Amey Pasarkar, and Itsik Pe'er Efficient and accurate inference of mixed microbial population trajectories from longitudinal count data. Cell Systems, 10(6):463-469 2020.
[37] J Atchison and Sheng M Shen. Logistic-normal distributions: Some properties and uses. Biometrika, 67(2):261-272 1980.
[38] Justin D Silverman, Heather K Durand, Rachael J Bloom, Sayan Mukherjee, and Lawrence A David, Dynamic linear models guide design and analysis of microbiota studies within artificial human guts. Microbiome, 6(1):1-20, 2018.
[39] Neal S Grantham, Yawen Guan, Brian J Reich, Elizabeth T Borer, and Kevin Gross. Mimix: bayesian mixed-effects model for microbiome data from designed experiments. Journal of the American Statistical Association, 115(530):599-60)) 2020.
[40] Dean Billheimer, Peter Guttorp, and William F Fagan. Statistical interpretation of species composition. Journal of the American statistical Association, 96(456):100-12214 2001.
[41] David M. Blei and John D. Lafferty. Correlated topic models. Advances Neural Information Processing Systems, pages 147-154, 2005.
[42] Vera Pawlowsky-Glan Glaahn, Juan José Egozcue, and Raimon Tolosana-Delgado. Modeling and analysis of compositional data. John Wiley & Sons, 2015.
[43] James T Morton, Jon Sanders, Robert Quinn, Daniel McDonald, Antonio Gonzalez, Yoshiki Vázquez-Baeza, Jose A Navas-Molina, Se Jin Song, Jessica L Metcalf, Embriette R Hyde, al. Balance trees reveal microbial niche differentiation. MSystems, 2(1), 2017.
[44] Justin D Silverman, Alex D Washburne, Sayan Mukherjee, and Lawrence A David. A phylo- genetic transform enhances analysis of compositional microbiota data. Elife, 6:e21887, 2017.
[45] Alex D Washburne, Justin Silverman, Jonathan W Leff, Dominic Bennett, John L Darcy, Sayan Mukherjee, Noah Fierer, and Lawrence A David. Phylogenetie factorization of compo- sitional data yields lineage-level associations in microbiome datasets. PeerJ, 5:e2969, 2017.
[46] Juan José Egozcue and Vera Sawllossk--Gahnn Groups of parts and their balances in compo- sitional data analysis. Mathematical Geology, 37(7):795-828 2005.
[47] George Cybenko. Approximation by superpositions of a sigmoidal function. Mathematics of control, signals and systems, 2(4):303-314 1989.
[48] Patrick Kidger and Terry Lyons. Universal approximation with deep narrow networks. In Conference on Learning Theory, pages 2306-2327, 2020.
[49] Stan modeling language users guide and reference manual.
[50] John C Gower. Generalized procrustes analysis. Psychometrika, 40(1):33-51, 1975.
[51] Danilo Jimenez Rezende, Shakir Mohamed, Ivo Danihelka, Karol Gregor, and Daan Wierstra. One-shot generalization in deep generative models. arXiv preprint arXiv:1603.05110 2016.
[52] Anupriya Tripathi, Alexey V Melnik, Jin Xue, Orit Poulsen, Michael J Meehan, Gregory Humphrey, Lingjing Jiang, Gail Ackermann, Daniel McDonald, Dan Zhou, et al. Intermit- tent hypoxia and hypercapnia, hallmark of obstructive sleep apnea, alters the gut microbiome and metabolome. MSystems, 3(3), 2018
[53] Shabnam Shalapour, Xue- Jia Lin, Ingmar N Bastian, John Brain, Alastair D Burt. Alexander A Aksenov, Alison F Vrbanac, Weihua Li, Andres Perkins, Takaji Matsutani, et al. Inflammation- induced iga+ cells dismantle anti-liver cancer immunity. Nature, 551(7680))340-345 2017.
12

---

Under review as a conference paper at ICLR 2021
[54] David Lovell, Vera Hewlowsky-Gllbbe Juan José Egozcue, Samuel Marguerat, and Jürg Bähler. Proportionality: a valid alternative to correlation for relative data. PLoS Comput Biol, 1033:e100075,,00
[55] Ionas Erb and Cedric Notredame. How should we measure proportionality on relative gene expression data? Theory in Biosciences, -35(1-2):1--3, 2016.
[56] Thomas P Quinn, Mark F Richardson, David Lovell, and Tamsyn M Crowley. propr: an r- package for identifying proportionally abundant features using compositional data analysis. Scientific reports, 7(1):1-9, 2017.
[57] John Aitchison, Carles Barceló- Vidal, José Antonio Martin-Feránndez, and Vera Pawlowsky- Glahn. Logratio analysis and compositional distance. Mathematical Geology, 32(3):27-2275, 2000.
[58] Cameron Martino, James Morton, Clarisse A Marotz, Luke Thompson, Anupriya Tripathi, Rob Knight, and Karsten Zengler. A novel sparse compositional technique reveals microbial perturbations. MSystems, 4(1), 2019.
[59] Cameron Martino, Liat Shenhav, Clarisse A Marotz, George Armstrong, Daniel McDonald, Yoshiki Vázquez-Baeza, James T Morton, Lingjing Jiang, Maria Gloria Dominguez-Ber Belll,, Austin D Swafford, al. Context-aware dimensionality reduction deconvolutes gut microbial community dynamics. Nature biotechnology, pages 1-4, 2020.
[60] Karl Pearson. Mathematical contributions to the theory of evolution.-on a form of spurious correlation which may arise when indices are used in the measurement of organs. Proceedings of the royal society of london, 60(359-367)))) (11-3366):44999999 1897.
[61] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al. Pytorch: An imperative style,
deep learning library. In Advances in neural information processing systems, pages 8026-8037) 2019.
[62] WA Falcon. Pytorch lightning. Gitllub. Note: http:///iihbbb.omm//yyyyothhhhnp lightning, 3, 2019.
[63] Daniel McDonald, Jose C Clemente, Justin Kuczynski, Jai Ram Rideout, Jesse Stombaugh, Doug Wendel, Andreas Wilke, Susan Huse, John Hufnagle, Folker Meyer, et al. The biological observation matrix (biom) format or: how i learned to stop worrying and love the ome-ome. Gigascience, 1(1:22047-217X, 2012.
[64] John D Hunter. Matplotlib: A 2d graphics environment. Computing in science & engineering, 93):90095 2007.
[65] Pauli Virtanen, Ralf Gommers, Travis E Oliphant, Matt Haberland, Tyler Reddy, David Cour- napeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, et al. Scipy 1.0: fundamental algorithms for scientific computing python. Nature methods, 17(3):261-722, 2020.
[66] Charles Harris, K Jarrod Millman, Stéfan J van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J Smith, et al. Array programming with numpy. arXiv preprint arXiv:2006. 10256, 2020.
[67] Max A Woodbury, Inverting modified matrices. Memorandum report, 42(106):336, 1950.
[68] William H Press, Saul A Teukolsky, Brian P Flannery, and William T Vetterling. Numerical recipes in Fortran 77: volume 1, volume of Fortrun numerical recipes: the art of scientific computing. Cambridge university press, 1992.
[69] Jonathan Chung-Kuse Huang and Tomasz Malisiewicz. Fitting a hierarchical logistic normal distribution. Unpublished Manuscript, pages 2-4.
13

---

Under review as a conference paper at ICLR 2021
[70] Antonio Gonzalez, Jose A Navas-Molina, Tomasz Kosciolek, Daniel McDonald, Yoshiki Vázquez-Bees -aazz, Gail Ackermann, Jeff DeReus, Stefan Janssen, Austin Swafford, Stephanie B Orchanian, et al. Qiit:: rapid, web-enabled microbiome meta-analysis. Nature methods, 15(10):7967798 2018.
[71] Daniel McDonald, Benjamin Kaehler, Antonio Gonzalez, Jeff DeReus, Gail Ackermann, Clarisse Marotz, Gavin Huttley, and Rob Knight. redbiom: a rapid sample discovery and feature characterization system. mSystems, 4(4), 2019.
[72] Torbjørn Rognes, Tomáš Flouri, Ben Nichols, Christopher Quince, and Frédéric Mahé. Vsearch: a versatile open source tool for metagenomics. PeerJ, 4:e2584, 2016.
Appendices
A CHALLENGES IN DERIVATION OF AN ANALYTICAL MULTINOMIAL VAE ELBO
Recall that the generative model for Multinomial PPCA is given as follows:
p(=1 v7) Mult((x/ (n))
(13)
p(n) 2) = (1) Wz u, (4,oIa-1)
(14)
p(=) = (z Ix)
(15)
With this in mind, we wish to estimate variational distributions q(n, x| = (< to ap- proximate the posterior p(n,, x)).. These variational distributions can both be chosen be normal distributions as follows:
(<=[>) = NVV (log(ar) 14), D)
<n)2) p(nl2)
Noting that NVVV' (log(z) x)) 11), D). ((]]))) can be derived from g(=[x) as follows:
<n]x)) = = WWV
μ1) + μ. DDW + +
To fine-tune these variational distributions to approximate the posterior distribution, we can mini- mize the following KL divergence:
arg max KL(g(n, =|=)||p((,,>>>>>>> 2|x))
<(n=(
Since cannot optimize this quantity directly, opt instead to maximize the evidence lower bound (ELBO) given by
1(7, 2|2)] ((1), 2|x)]
2 Eg( log p(zl p(z(m))(()a)pp())
Bg(n))))))
log
q 217))
We can partition this lower bound into three parts, given as follows:
p(n|=))
[x(x(0) log p(z)
Ba Dog (n(()) log
((nl)))
(61
(116)
Since Mult(æ/p) x Σ x. log Ps. the first term is given by
[logp(>|)) (10ola)gsste((00p((a111) & Ex(()
Σ Ro Un)i
f N NWW (log(æ) 1)), WDW + +021) E (nn) dn
14

---

Under review as a conference paper at ICLR 2021
Estimating the above integral is equivalent estimating the expectation of logistic normal distri- bution, which does not have an analytical solution (37) As a result, the analytical ELBO for the Multinomial VAE is intractable.
B THE ILR TRANSFORM
As outlined in Equation any orthogonal basis T can be used to perform the ILR transform. However, there are select few bases that can represented by binary tree. To see how an orthogonal basis can be constructed from a binary tree, consider the illustration in Figure
771
1/2
1)3
+1
-1
-1
-1
7/1
+1
+1
-1
7/2
+1
-1
0
1/3
X1
X2
X3
X44
Figure S1: A small example showing how an orthogonal basis can be constructed from binary tree. = (I1, 22, 33, 33) E S¹ represents species proportions and = (m 72, 7)s) E R represents the log-ratios in the internal nodes. V matrix of orthogonal contrasts that also can be represented from binary tree.
Here the rows of the matrix in Figure $1 orthogonal and the resulting product v T log yields
F243
T3
log
1/2
log
T/3
log
274
.74
The contrast matrix can forced be orthonormal such that & T
Ia-1. 1. as highlighted in
Equation 6 Furthermore, this construction can be scaled to large binary trees as shown in Figure $2 Here, 7) represents the log-ratios at the internal node given by
|
gxxr) log
7
(16)
g (xx
The runtime of this operation is further discussed in Appendix B.I
B. THE RUNTIME OF THE ILR TRANSFORM
The naive runtime of the ILR transform of single sample E Rd ((2) due to the running time of dense matrix-vector multiplication. As shown in Equation 5 the ILR transform can also be represented by log-linear transformation with a contrast matrix. The binary tree can used to represent contrast matrix, discussed in (46).
the binary tree balanced, each row of Y will have O(log d) non-zero elements, since the tree has height of Ollog d). Given that there are rows, the matrix- multiplication behind the inverse ILR transform in can be done in (d log d). For the same reason, the ILR transform has a runtime of Odlogd).
15

---

Under review as a conference paper at ICLR 2021
n
k
r
S
t
Figure S2: An illustration of how the ILR basis can be constructed on large trees. The quantities g(x,)) and g(x.) yield the geometric means within a vector of proportion for subsets æ, and æ Here, r and refer to the sets of features in the left and right subtrees for the internal node I. The log-ration 77 can be obtained from either Equation or Equation [6
B.2 SHIFT INVARIANCE OF THE SOFTMAX TRANSFORM
There are two different scale identifiability issues. The softmax transform that is commonly used shift invariant, where
<(x) + oxx Va E R
This shift invariance will cause an identifiability issue when identifying the decoder matrix W. The inverse ALR transform resolves this issue by enforcing an isomorphism; of the coordinates is set to zero as follows:
ALR(æ)= = log I1
log
(x) ((11....,,5 FT-1,,
Jrd
Fd
One of the implicit constructions of this transform that the resulting contrast matrix is not or- thogonal (42). Furthermore, because the ALR transform is isometric, the resulting Euclidean distances in will not approximate the Aitchison distance. The ILR transform provides the best of both worlds, enforeing both isometry and isomorphism. Furthermore, there is a tight connection between the ILR transform and compositional PCA that is further discussed in the next section.
B.3 GLOBAL LOG-CONVEXITY OF THE MULTINOMIAL LOGISTIC NORMAL DISTRIBUTION
Since the multinomial logistic normal distribution difficult to directly evaluate, challenging to make statements regarding its maximum likelihood estimator. With the posterior factorization highlighted Equation 4 we can make concrete statements about the original posterior factors. The log probability density of the multinomial distribution p(z|)) can be written
logp(a))) x
T; log di (n))
F; (T)); - m log (n))
(()) T(x:) A(n)
where m Σ It is the total number of counts and the functions g(r);: (Mr)) T(x); Zi. and A(n) = log E/= exp((X)) are the natural parameters of the exponential family distribu-
16

---

Under review as a conference paper at ICLR 2021
tion. The Hessian of log p(z)n) is given by d logp(zln) dnjdnn
d2 dAA(n)
Since A(n) is strietly convex, log (7) is also strictly convex.
Similarly, the log probability density of the multivariate Gaussian distribution p(z|n) = W also strictly convex with respect to f and E. Since the multinomial logistic normal distribution can be written the sum of log p(=| (==))) m) and log p(alm)) it is also strictly convex. Therefore, there must be a unique optimal estimate for μ and E.
B.4 MULTINOMIAL VARIATIONAL AUTOENCODER ESTIMATION
The covariance matrix W WDWT + I E Rd--xxd1 in Equation 11 can be efficiently inverted using the Woodbury identity (67 68) The prior (z is evaluated given log (z) = ( 00 Ix), where ( = Fr. T (log( x (1))
H)) obtained through gradient descent optimiza- tion. Once (1) is obtained, the remaining VAE parameters 0 can be estimated via gradient descent optimization. Like the EM algorithm proposed in (69), this procedure alternates between estimat- ing (1) and O and repeats until convergence. For single-layer linear encoder, this optimization procedure will eventually reach the global maxima with respect to H and the multivariate normal mean and covariance in g( H 10, x()), due to the log-convexity of the multinomial and normal distributions (Appendix B.3 On fully observed data, W and D can be estimated up to rotation and scale (Appendix B.6). This is not guaranteed for sparse count data, but using simulations, we can empirically show that increasing the complexity encoder architectures can help. Since H needs be accurately estimated before optimizing Fz, W, a, and D, multiple gradient descent steps are required for given mini-batch.
The optimal (t) that maximizes log ((x (((( (1)
In practice obtaining the optimal (1) for a given mini-batch i may require hundreds of gradient descent steps from random initialization. Since (i) is used to approximate the multinomial parameters describing the observed counts æ, initialize H with ( ' T X ()). In practice, this can greatly reduce the number of gradient descent updates needed per mini-batch.
B.5 THE CONNECTION BETWEEN MULTINOMIAL VAEs AND COMPOSITIONAL PCA
With singular value decomposition, the resulting factors can be used to approximate the row and column distances. For singular value decomposition given by X USVT, distances and column distances can be approximated as follows:
[x. æj. l Rs || - 8ytt;11
| x e.ile |s,0 8y0;112
This relationship also yields connection the row and column covariances; the row covariances are given by X U and the column covariances are given by XX" VS
A similar relationship applies to compositional PCA, except the singular value decomposition is applied to CLR-transformed values of X. The CLR transform is given by
CLR(#) = log F1
log
Zd
CLR (x) 6(x)
g(x)
The inverse CLR transform is equivalent to the softmax transform Φ(x). The Euclidean distance on CLR and ILR transformed values is given by the Aitchison distance (57). As result, the row and column distances are given by the Aitchison distance. The covariance matrix of CLR- RuR-trranffrmedd values is singular, and as a result, the resulting singular value decomposition will have full rank of d-1. As shown by Egozcue et al. (23), the top d- components of the singular value decomposition are all in ILR coordinates. This observation cements the connection between compositional PCA and our proposed Multinomial VAE; since all of the singular value components be represented
17

---

Under review as a conference paper at ICLR 2021
in ILR coordinates, this theoretically justifies the use of the ILR transform within proposed Multinomial VAE.
Due to this connection, the distances between the learned VAE representations across pairs of sam- ples should be proportional to the Aitchison distance. Furthermore, the distances between rows of W would also be given by the Aitchison distance, which is equivalent to Lovell's mn-nommalized proportionality constant (23).
B.6 LACK OF SCALE IDENTIFIABIITT OF PPCA
As discussed in (8) VAEs will have the same identifiaiility issues that PPCA has; namely, for a diagonal matrix A Biok the following equivalences hold:
W WA. V+
As result, the decoder weights W can only be identified up to scale. Furthermore, the Multinomial VAE architecture is asymmetric and, as result, the Tranpose Theorem proposed by (17) not as readily applicable. The lack of symmetry between W and V complicates the process of obtaining cigenvalues from W.
C SIMULATION AND MICROBIOME BENCHMARKS
C.I SIMULATION DETAILS
In the dense simulation highlighted Figure there were 100 features, 1000 samples and million counts per sample. These counts were drawn from multinomial logistic normal distribution whose covariance matrix had a rank of 10. The sparse dataset highlighted dimensions and 1000 samples, each of which contained 100 counts that were sampled from a multi- nomial logistic normal distribution whose covariance matrix had a rank 10. The sparse dataset highlighted in Figure -p contained 5000 dimensions and 10000 samples, each of which contained only 1000 counts total. These counts were sampled from a multinomial logistic normal distribution whose covariance matrix had rank of 128.
Figure 2a-h contained 200
The encoder architecture used for the stochastic and analytical VAEs consisted of fully connected dense layers with interwoven softplus activation functions. The intermediate encoder layers all have the same input and output dimension as the number of latent dimensions specified in the model. The same architecture was used for both the stochastic VAEs and the analytical VAEs. A schematic of the architecture is shown in Figure S3
Multinomial PCA was fitted using both penalized likelihood estimation and Hamiltonian Monte Carlo using Stan (Appendix (C.).. The initialization of these models was determined by the es- timated factor loadings from the analytical VAE. This particularly critical for fitting HMC on high-dimensional datasets. We were not able to run Stan on the sparse datasets.
C.2 EXPERIMENTAL DATASET DETAILS
The full 8.8k mouse pellet dataset was retrieved from qïita (70) using redbiom (7I)
In the hypoxia study conducted by Tripathi et al., mice were placed simulated conditions to induce intermittent hypoxic/hypercapnic (IHH) stress, where the mouse's oxygen supply was reduced and its CO2 supply was increased. The goal of this experiment was detect the gut microbiome differ- ences between the IHH mice and the control mice. We will use this information as a classification benchmark to determine the experimental conditions the mice were placed under can be predicted from 16S sequencing counts. In total, there were 48 mice: 24 of these mice were placed in these IHH conditions and the remaining 24 mice served as controls. Fecal samples were collected twice a week for 6 weeks, resulting in 579 samples total. There were 5775 microbial taxa that were detected in this dataset.
The goal of the hepatocellular carcinoma (HCC) study conducted by Shalapour was to investigate the interplay between immunity, diet and carcinoma. In total, there were 478 mice with 52 phenotypes, split according to diet, immunity status and whether not the mice were induced with HCC. Our
18

---

Under review as a conference paper at ICLR 2021
Encoder
Decoder
x
7
X
q (z
p(llz) (xxln)
Figure S3: Visualization of VAE architecture
classification objective to predict the mouse phenotype from the 16S sequencing counts obtained from their fecal pellets. In total, there was one fecal pellet collected from each mouse and there were 2794 mierobial taxa detected across these samples. This benchmark is more challenging than the IHH benchmark due to the extreme class imbalance; some phenotypes contain upwards of 30 mice and other contain as few as mice (Figure S4). The resulting classification benchmarks are shown in Table
C.3 MICROBIOME VAE TRAINING DETAILS
For the full 11k fecal pellet dataset, we focused on samples that were processed using closed- reference OTU picking using vsearch (72). This dataset was split into 80/10/10 cain/validation/test split in order measure out-of-distribution generalizability.
We trained our proposed analytical VAE and the stochastic VAE with 128 latent dimensions and batch size 1000 samples for 10000 epochs, with a learning rate of 10 and 100 gradient descent steps per mini-bateh Checkpoints were recorded every epoch; the models with the best validation error are reported. Cosine annealing with warm restarts was used as a learning rate scheduler, with the intention of easily escaping saddle points during optimization.
C.4 METRICS FOR EVALUTING ESTIMATED DECODER WEIGHTS
To measure the agreement between the ground truth simulations and the estimated decoder weights W, three different metrics evaluated, namely, axis alignment and subspace distance, as proposed in (18), in addition to Procrustes analysis (50). Given the decoder weights W and the ground truth principal components U, these metrics are defined follows:
Axis alignment
1 h
(UWW)) )
d(U. w)
max
U.
WiIlE
This metric is ultimately a measure the average cosine distance between the ground truth eigen- vectors and the estimated decoder weights.
19

---

Under review as a conference paper at ICLR 2021
Subspace distance
W) -
where W. denotes the left singular vectors of W, given by W W.AV. Since U yields the ground truth correlations and W. W. yields the estimated correlations, this metric can be interpreted as measure of agreement between the ground truth correlations and the estimated correlations.
Procrustes Analysis
d(U, W) = arggginn m W WRAll R,A
for A, R E R where A is diagonal matrix and R is rotation matrix. Prior evaluating this metric, both U and W are standardized such that Tr(U'²)) = Tr(WWT) = and W is centered around the origin.
Correlation
d(, w) = Corr(vee(Au), vec(Aw))
where (ir = I|u ullll and Aw (i.j) = I||w; - wy /Il2 denote the pairwise distances of U and W. The pairwise distances rotation invariant. Furthermore, the measure of the correlation is agnostic to scale, and will ignore the eigenvalue scale identifiability issues highlighted in the main text.
C.5 STAN MULTINOMIAL PCA IMPLEMENTATION
data {
Int<lower=0> N;
1/ number of samples
int<lower=0> D;
// number of dimensions
int<lloer=0>> K;
// number of latent dimensions
matrix (D-1, D] Psi;
// Orthonormal basis
int y(N, D]:
/1 observed counts
parameters matrix [N, D-1] eta; // ilr transformed abundances
matrix [D-1, K] W;
real <law<roo sigma;
transformed parameters { matrix [D-1, 1, D-1] 1] Sigma;
matrix [D- 1, D-1] I;
vector[D [D-1] z:
[lag_matrix (rep_veetor (1.0, D)))); Sigma = W W + square (sigma) I;
z = rep-vector (0, D-1);
model {
// generating counts
for (n in 1:NH
eta (ecieniallSoofffofa wowftmaxtto.vttaaooooooooootetee (to_vector [n] Psi )))
(z, Sigma);
y[n] multinomial
}
20

---

Under review as a conference paper at ICLR 2021
HCC Phenotype Distribution
Number of Mice
30
20
10
0
Phenotype
Figure S4: Distribution of 52 phenotypes the Shalapour study, highlighting the class imbal- ance inherent in this benchmark.
Dataset
Method corrected
Pred Log prob
score
AUC
Shalapour et al.
analytical- VAE
57750
0.845
0.820
stochastic-Val -VA
11615
1.000
1.000
LDA raw
169384
0.800
0.988
NA
1.000
1.000
Tripathi et
analytical- VAE
117789
0.845
0.820
(tochastic-AA
13815
0.914
0.901
LDA raw
132091
0.667
0.671
NA
0.571
0.632
All Test Samples
analytical-VAr stochastic- VAE
117789
13647
LDA
169384
Table 2: Classification on held out data benchmarks. KNN classifiers were trained on data observed by the VAE, and evaluated on 1149 held out test samples. KNN with k=1 was used across all representations. Predictive log probability is also added to measure the generalizability of the VAEs across held out samples (smaller values are better).
Due to the relationship between squared Euclidean distance and proportionality given in Equation 12 we expect the proportionality metric to exhibit log-linear relationship between the Euclidean distances obtained from the VAE embeddings. Indeed, this is apparent with the analytical VAE embeddings in Shalapour et al., and to lesser extent in Tripathi et al., as shown Figure S5
21

---

Under review as a conference paper at ICLR 2021
Shalapour et al.
Log Proportionality
3
3
2
2
1
0
-1
-1
-2
-1
0
3
4
5
Log Embed Distance
Log Embed distance
Tripathi et al.
Log Proportionality
2
2
0
-2
-2
-1
0
2
4
Log Embed Distance
Log Embed distance
Figure S5: Comparison of Lovell's proportionality metric and the pairwise embedding distances for the stochastic VAE and the analytical VAE evaluated on the Shalapour al. dataset and the Tripathi et al. dataset. The embedding distance is defined by the Euclidean distance between two rows within WW Since the proportionality metric is not defined for zeros, it is only evaluated for samples where both taxa are observed. Only the top 300 most abundant microbes visualized here.
22


Reviews:
[{"review_id": "T6QXUWxZZLH", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper extends prior results, namely that VAEs are able to learn the principal components. The novelty is the extension to a new distribution: multinomial logistic-normal distribution. This is achieved by using the Isometric log-ratio (ILR) transform. While prior results were derived analytically, this paper provides (only) empirical evidence for the claim regarding the multinomial logistic-normal distribution. \n\nOverall, I don’t find that the provided experiments provide convincing support of  the claim of the paper, namely that VAEs are able to learn the principle components for the multinomial logistic-normal distribution.\n\nWhile the proposed approach yields better results than alternative approaches in Figure 1, it is not clear to me why the shown results indicate that the VAE was actually able to learn the principal components: based on my understanding of the used metrics (for instance,  axis-alignment of 0 indicates perfect results, and 1 is the worst): the proposed approach achieves 0.8 on 200 dimension and 0.95 on 1000 dimensions. When a deeper model is used in Figure 1, the metrics go down, which is good, but they still stay above 0.5, i.e., far away from 0. Based on my understanding of the used metrics, I am not sure why this would provide empirical support of  the claim in this paper that the principal components can be recovered. \n\nIn the definition of axis alignment in the appendix, it seems like the enumerator is missing a square, compared to the original definition in the referenced paper [14].\n\nThe authors also claim that Figure 2 shows that the learned embedding dimensions are orthogonal to each other. While I can believe this in the right plot, I am not convinced by the left plot, where a different scaling is used (as mentioned by the authors), and hence the diagonal does not seem so much larger than the off-diagonals to me.\n\nThe paper also proposes a batch-correction as additional improvement, but Table 1 shows that it makes results  actually worse. \n\nApart from that, Table 1 would also benefit from adding some state-of-the-art baseline models, e.g., as discussed in the related works section. Simply applying KNN to the raw counts seems like a very basic baseline.\n\nApart from that, the writing of the paper could be improved a lot.  I found it very confusing to come across references to Figure 5 etc,  and then not finding them in the paper ... eventually I realized that they were in a separate supplement. This could be made clearer in the paper. Moreover, there are also several typos in the paper (e.g., “thethe”), and some grammatically incorrect sentences. Also,  in the related works section, the authors refer to work by  “Luca et al” without a citation--does it  refer to Lucas (with an s) [4] ?\n\t\t\n++++ updates after discussion period ++++\n\nWhile the initial paper provided empirical evidence  for the claims of the paper, based on the reviews an appendix of about 9 pages was added in the revised version, and the focus now seems on providing theoretical support of the claims in the paper. Also several experimental results were added in the main paper in response to the reviews. I feel like all these (quite major) changes indicate that this paper is not mature enough for publication at this point. So I maintain my current review.", "extra_scores": {}, "extra_text": {}}, {"review_id": "FcuvI6zUkY7", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper is well written, presenting a great interdisciplinary work on covariance estimation. It relies on recent techniques such as connection between VAE and PPCA, ILR transformation, and presents an augmented VAE to obtain MAP estimates on multinomially distributed data. Such setting has direct application potential in many bioinformatics problems. The methodology is quite dense which reads look to me (though not in every detail). \n\nThe only chance it gets rejected is its relatively narrow audience group in ICLR. \n\n1. the notation can be improved to fit community convention: say P(z|\\theta, x), it's better to read P(z|x; \\theta) since \\theta are only parameters, not the target random variable.\n2. Clarity can be improved on the part elaborating Algorithm 1, to clarify conclusive statements, eg justification on complexity, and put aside technical details. Readers may have confusions on how ILR transformation is used here: maybe Section 4.2 can hint on the exact use case and the part after Algorithm 1 can be better organized.\n3. Experiments part should have more clear intro on machine learning abstracted versions of the problems, which better fit general audience in ICLR.\n4. Some minor typos, eg Section 4.1 Equation (2) should be I_{d-1} instead of I_d.\n\nHope this paper gets attention and wide adoptions on related biological applications.", "extra_scores": {}, "extra_text": {}}, {"review_id": "CVLdQJCHIc", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "### Summary\nThe authors propose a framework to estimate high-dimensional covariance matrices using Multinomial variational autoencoders, showing application to biomedical data sets.\nIn particular, they use the probabilistic PCA framework, extending it to multinomial-distributed data. They show that similar work on PPCA are limited to Gaussian data, and non-easily applicable to bag-of-words. Analysis on real-world biomedical data sets are promising.\n\n### Reasons for score\nI fail to grasps the novelty with respect to the state of the art. The authors cite numerous of similar works, such as mixture modelling (LDA). The authors mention that such techniques rely on Dirichlet distributions, but that is not always the case (such as in [1] and [2]), where they show it can be approximated with a Logistic Normal. It seems such models can be used in the applications mentioned by the authors. However, since there are no comparison in the paper, I fail to understand the benefit of one with respect to the others.\n\n### Pros\n- ILR transform to deal with identifiability of softmax\n- Analysis on synthetic and real data sets\n\n### Cons\n- Technical details can be improved. For example, it is not straightforward to understand what the encoder / decoder parameters refer to (is the type of activation functions at each layer? dimension of layers? number of layers?) \n- There is no comparison with mixture models that infer covariance matrices (like [1] and [2]). It is difficult to assess what are the benefit of one with respect to the others.\n\n### Minors\n- \"thethe\", page 5\n\n### References\n[1] Srivastava, Akash, and Charles Sutton. \"Autoencoding variational inference for topic models.\" arXiv preprint arXiv:1703.01488 (2017).\n[2] Blei, David, and John Lafferty. \"Correlated topic models.\" Advances in neural information processing systems 18 (2006): 147.", "extra_scores": {}, "extra_text": {}}, {"review_id": "8TbGG_y0451", "reviewer": "AnonReviewer4", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The authors demonstrate a VAE model and estimation framework with which the PPCA subspace is recovered for data with multinomial observations. The authors specifically are interested in the high-d scenario and show that their analytical results out-perform other VAE methods. \n\n## General comments\nThe paper offers a novel solution to the seemingly persistent problem of performing efficient inference of a non-congugate latent variable model for count data. The authors claim that existing solutions are not scalable. The paper is clearly written with a few minor typos. I have 2 complaints: \n1) The benchmarking analysis reports only relative performance across other VAEs. There are no other methods even considered and it would be nice to have some sense of what a \"good\" result looks like in an absolute sense (eg. what's a good score on subspace distance for 200 or 1000 input dimension problem? What would you get if you didn't use the VAE approach?).  \n2) Could the authors comment more on the scalability problems they are talking about? There appear to be a number of scalable PCA techniques available where the full covariance matrix need not be formed and decomposed in memory. Would none of those approaches  be adaptable the their problem?\n\n## minor note\nThe authors mention a few papers where MCMC was proposed as a means to obtain a posterior over the latent variables but these ventures were dismissed as suffering from the curse of dimensionality (or something to the effect of that argument). One outstanding oversight to this collection of papers is [Linderman et al.](https://papers.nips.cc/paper/5660-dependent-multinomial-models-made-easy-stick-breaking-with-the-polya-gamma-augmentation) wherein the stick-breaking construction of the multinomial is coupled with Polya-gamma augmentation to offer a highly efficient sampling procedure for precisely the model the authors are proposing. Could the authors comment in the rebuttal on how this approach ", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: Authors extend the probabilistic PCA framework to multinomial-distributed data. Scalable estimation of principal components in the model is achieved using a multinomial variational autoencoder in combination with an isometric log-ratio (ILR) transform.
The reviewers did not agree on the degree of novelty of the paper to PC estimation.
The presentation of the paper can be improved.
The reviewers criticise that large changes have been made to the paper during the rebuttal phase.
Overall, the paper is borderline and due to the mentioned large changes I recommend a rejection (and re-review at a different venue).

Decision: Reject