Year: 2020

Abstract:
Variational autoencoder (VAE) estimates the posterior parameters (mean and variance) of latent variables corresponding to each input data. While it is used for many tasks, the transparency of the model is still an underlying issue. This paper provides a quantitative understanding of VAE property by interpreting VAE as a non-linearly scaled isometric embedding. According to the Rate-distortion theory, the optimal transform coding is achieved by using a PCA-like orthonormal transform where the transform space is isometric to the input. From this analogy, we show theoretically and experimentally that VAE can be mapped to an implicit isometric embedding with a scale factor derived from the posterior parameter. As a result, we can estimate the data probabilities in the input space from the prior, loss metrics, and corresponding posterior parameters. In addition, the quantitative importance of each latent variable can be evaluated like the eigenvalue of PCA.

heading: 1 INTRODUCTION
text: Variational autoencoder (VAE) (Kingma & Welling, 2014) is one of the most successful generative models, estimating posterior parameters of latent variables for each input data. In VAE, the latent representation is obtained by maximizing an evidence lower bound (ELBO). A number of studies (Higgins et al., 2017; Kim & Mnih, 2018; Lopez et al., 2018; Chen et al., 2018; Locatello et al., 2019; Alemi et al., 2018; Rolı́nek et al., 2019) have tried to reveal the property of latent variables. However, quantitative behavior of VAE is still not well clarified. For example, there has not been a theoretical formulation of the reconstruction loss and KL divergence in ELBO after optimization. More specifically, although the conditional distribution pθ(x|z) in the reconstruction loss of ELBO is predetermined such as the Gaussian or Bernoulli distributions, it has not been discussed well whether the true conditional distribution after optimization matches the predetermined distribution.
Rate-distortion (RD) theory (Berger, 1971), which is an important part of Shannon information theory and successfully applied to image compression, quantitatively formulates the RD trade-off optimum in lossy compression. To realize a quantitative data analysis, Rate-distortion (RD) theory based autoencoder, RaDOGAGA (Kato et al., 2020), has been proposed with isometric embedding (Han & Hong, 2006) where the distance between arbitrary two points of input space in a given metrics is always the same as L2 distance in the isometric embedding space. In this paper, by mapping VAE latent space to an implicit isometric space like RaDOGAGA on variable-by-variable basis and analysing VAE quantitatively as a well-examined lossy compression, we thoroughly clarify the quantitative properties of VAE theoretically and experimentally as follows.
1) Implicit isometric embedding is derived in the loss metric defined space such that the entropy of data representation becomes minimum. A scaling factor between the VAE latent space and implicit isometric space is formulated by the posterior for each input. In the case of β-VAE, the posterior variance of each dimensional component in the implicit isometric embedding space is a constant β/2, which is analogous to the rate-distortion optimal of transform coding in RD theory. As a result, the reconstruction loss and KL divergence in ELBO can be quantitatively formulated. 2) From these properties, VAE can provide a practical quantitative analysis of input data. First, the data probabilities in the input space can be estimated from the prior, loss metric, and posterior parameters. In addition, the quantitative importance of each latent variable, analogous to the eigenvalue of PCA, can be evaluated from the posterior variance of VAE.
This work will lead the information theoretic generative models in the right direction.

heading: 2 RELATED WORKS

heading: 2.1 VARIATIONAL AUTOENCODER AND THEORETICAL ANALYSIS
text: In VAE, ELBO is maximized instead of maximizing the log-likelihood directly. Let x ∈ Rm be a point in a dataset. The original VAE model consists of a latent variable with fixed prior z ∼ p(z) = N (z; 0, In) ∈ Rn, a parametric encoder Encφ : x ⇒ z, and a parametric decoder Decθ : z ⇒ x̂. In the encoder, qφ(z|x) = N (z;µ(x),σ(x)) is provided by estimating parameters µ(x) and σ(x). Let Lx be a local cost at data x. Then, ELBO is described by
ELBO = Ex∼p(x) [ Ez∼qφ(z|x)[log pθ(x|z)]−DKL(qφ(z|x)‖p(z)) ] . (1)
In Ex∼p(x)[ · ], the first term Ez∼qφ(z|x)[ · ] is called the reconstruction loss. The second term DKL(·) is a Kullback–Leibler (KL) divergence. Let µj(x), σj(x), and DKLj(x) be j-th dimensional values of µ(x), σ(x), and KL divergence. Then DKL(·) is derived as:
DKL(·) = n∑ j=1 DKLj(x), where DKLj(x) = 1 2 ( µj(x) 2 + σj(x) 2 − log σj(x)2 − 1 ) . (2)
D(x, x̂) denotes a metric such as sum square error (SSE) and binary cross-entropy (BCE) as loglikelihoods of Gaussian and Bernoulli distributions, respectively. In training VAE, the next objective is used instead of Eq. 1, where β is a parameter to control the trade-off (Higgins et al., 2017).
Lx = Ez∼qφ(z|x)[D(x, x̂)] + βDKL(·). (3) However, it has not been fully discussed whether the true conditional distribution matches the predetermined distribution, or how the value of KL divergence is derived after training.
There have been several studies to analyse VAE theoretically. Alemi et al. (2018) introduced the RD trade-off based on the information-theoretic framework to analyse β-VAE. However, they did not clarify the quantitative property after optimization. Dai et al. (2018) showed that VAE restricted as a linear transform can be considered as a robust PCA. However, their model has a limitation for the analysis on each latent variable basis because of the linearity assumption. Rolı́nek et al. (2019) showed that the Jacobian matrix of VAE at each latent variable is orthogonal, which makes latent variables disentangled implicitly. However, they do not uncover the orthonormality and quantitative properties because they simplify KL divergence as a constant. Dai & Wipf (2019) also showed that the expected rate of VAE for the r-dimensional manifold is close to −(r/2) log γ +O(1) at γ → 0 when pθ(x̂|x) = N (x̂;x, γIm) holds. The remaining challenge is to clearly figure out what latent space is obtained at a given dataset, a loss metric, and β in the model.

heading: 2.2 RATE-DISTORTION THEORY, TRANSFORM CODING, AND ISOMETRIC EMBEDDING
text: RD theory (Berger, 1971) formulated the optimal transform coding (Goyal, 2001) for the Gaussian source with square error metric as follows. Let x ∈ Rm be a point in a dataset. First, the data are transformed deterministically with the orthonormal transform (orthogonal and unit norm) such as Karhunen-Loève transform (KLT) (Rao & Yip, 2000). Let z ∈ Rm be a point transformed from x. Then, z is entropy-coded by allowing equivalent stochastic distortion (or posterior with constant variance) in each dimension. A lower bound of a rate R at a distortion D is denoted by R(D). The derivation of R(D) is as follows. Let zj be the j-th dimensional component of z and σzj2 be the variance of zj in a dataset. It is noted that σzj2 is the equivalent to eigenvalues of PCA for the dataset. Let d be a distortion equally allowed in each dimensional channel. At the optimal condition, the distortion Dopt and rate Ropt on the curve R(D) is calculated as a function of d:
Ropt = 1
2 m∑ j=1 max(log(σzj 2/d), 0), Dopt = m∑ j=1 min(d, σzj 2). (4)
The simplest way to allow equivalent distortion is to use a uniform quantization (Goyal, 2001). Let T be a quantization step, and round(·) be a round function. Quantized value ẑj is derived as kT , where k = round(zj/T ). Then, d is approximated by T 2/12 as explained in Appendix H.1.
To practically achieve the best RD trade-off in image compression, rate-distortion optimization (RDO) has also been widely used (Sullivan & Wiegand, 1998). In RDO, the best trade-off is
achieved by finding a encoding parameter that minimizes a cost L = D + λR at given Lagrange parameter λ. Recently, deep image compression (Ballé et al., 2018) has been proposed. In these works, instead of an orthonormal transform with sum square error (SSE) metric in the conventional lossy compression, a deep autoencoder is trained with flexible metrics, such as structural similarity (SSIM) (Wang et al., 2001) for RDO. Recently, an isometric autoencoder, RaDOGAGA (Kato et al., 2020) was proposed based on Ballé et al. (2018). They proved that the latent space to be isometric to the input space if the model is trained by RDO using a parametric prior and posterior with constant variance. By contrast, VAE uses a fixed prior with a variable posterior. In section 3, we explain that VAE can be quantitatively understood as the rate-distortion optimum as in Eq. 4 by mapping VAE latent space to implicit isometric embedding on a variable-to-variable basis as in Fig. 1 .

heading: 3 UNDERSTANDING OF VAE AS A SCALED ISOMETRIC EMBEDDING
text: This section shows the quantitative understanding of VAE. First, we present the hypothesis of mapping VAE latent space to an implicit isometric embedding space. Second, we reformulate the objective of β-VAE for easy analysis. Third, we prove the hypothesis from the minimum condition of the objective. Then, we show that ELBO can be interpreted as an optimized RDO cost of transform coding where the quantitative properties are well clarified, as well as discuss and correct some prior theoretical studies. Lastly, we explain the quantitative properties of VAE to validate the theory including approximations and provide a practical data analysis.

heading: 3.1 HYPOTHESIS OF MAPPING VAE TO THE IMPLICIT ORTHONORMAL TRANSFORM
text: Figure 1 shows the mapping of VAE to the implicit isometric embedding. Assume the data manifold is smooth and differentiable. Let Sinput(⊂ Rm) be an input space of the dataset. D(x, x́) denotes a metric for points x, x́ ∈ Sinput. Using the second order Taylor expansion, D(x,x + δx) can be approximated by tδx Gxδx, where Gx and δx are an x dependent positive definite Hermitian metric tensor and an arbitrary infinitesimal displacement in Sinput, respectively. The derivations of Gx for SSE, BCE, and SSIM are shown in Appendix H.2. Next, an implicit isometric embedding space SIso(⊂ Rm) is introduced like the isometric latent space in RaDOGAGA (Kato et al., 2020), such that the entropy of data representation is minimum in the inner product space ofGx. Let y and yj be a point in SIso and its j-th component, respectively. Because of the isometricity, p(x) ' p(y) will hold. We will also show the posterior variance of each dimensional component yj is a constant β/2. In addition, the variance of yj will show the importance like PCA when the data manifold has a disentangled feature by nature in the metric space ofGx and the prior covariance is diagonal.
Then, SIso is nonlinearly scaled to the VAE’s anisometric orthogonal space SVAE(⊂ Rn) on a variable-by-variable basis. Let z be a point in SVAE, and zj denotes the j-th component of z. Let p(yj) and p(zj) be the probability distribution of the j-th variable in SIso and SVAE. Each variable yj is nonlinearly scaled to zj , such that dzj/dyj = p(yj)/p(zj) to fit the cumulative distribution. dzj/dyj is σj(x)/ √ β/2, the ratio of posterior’s standard deviations for zj and yj , such that KL divergences in both spaces are equivalent. In addition, dimensional components whose KL divergences are zero can be discarded because such dimensions have no information.
3.2 REFORMULATION OF OBJECTIVE TO THE FORM USING ∂x/∂zj AND ∂x/∂z
We reformulate the objective Lx to the form using ∂x/∂zj and ∂x/∂z. Here, the dimensions of x and z, i.e., m and n, are set as the same. The condition to reduce n is shown in section 3.3.
Reformulation of D(x, x̂) loss: In accordance with Kato et al. (2020), the loss D(x, x̂) can be decomposed into D(x̆, x̂) + D(x, x̆), where x̆ denotes Decθ(µ(x)). The first term D(x̆, x̂) is a distortion between the decoded values of µ(x) with and without noise σ(x). We call this term as a coding loss. This term is expanded as follows. δx̆ denotes x̂ − x̆. Then, D(x̆, x̂) term can be approximated by tδx̆ Gxδx̆. Let xzj be ∂x/∂zj at zj = µj(x), and δzj ∼ N (0, σj(x)) be an added noise in zj . Then, δx̆ is approximated by δx̆ ' ∑m j=1 δzj xzj . Because δzj and δzk for j 6= k are uncorrelated, the average of D(x̆, x̂) over z ∼ qφ(z|x) can be finally reformulated by
Ez∼qφ(z|x) [D(x̆, x̂)] ' Ez∼qφ(z|x) [ tδx̆ Gxδx̆ ] ' n∑ j=1 σj(x) 2 txzjGxxzj . (5)
The second term D(x, x̆) is a loss between the input data and Decθ(µ(x)). We call this term a transform loss. We presume VAE is analogous to the Wiener filter (Wiener, 1964; Jin et al., 2003) where the coding loss is regarded as an added noise. From the Wiener filter theory, the ratio between the transform loss and coding loss is close to the ratio between the coding loss and the variance of the input data. The coding loss, approximately nβ/2 as in Eq. 14, should be smaller than the variance of the input data to capture meaningful information. Thus the transform loss, usually small, is not considered in the following discussion. Appendix B explains the detail in a simple 1-dimensional VAE. We show the exhaustive and quantitative evaluation of coding loss and transform loss in the toy dataset in appendix E.2 to validate this approximation.
Reformulation of KL divergence: When σj(x) 1, σj(x)2 − log σj(x)2 is observed. For example, when σj(x)2 < 0.1, we have −(σj(x)2/ log σj(x)2) < 0.05. In such dimensions, DKLj(x) can be approximated as Eq. 6 by ignoring the σj(x)2 term and setting p(µj(x)) to N (zj ; 0, 1):
DKLj(x) ' 1
2
( µj(x) 2 − log σj(x)2 − 1 ) = − log ( σj(x) p(µj(x)) ) − log 2πe
2 . (6)
Eq. 6 can be considered as a rate of entropy coding for a symbol with mean µj(x) allowing quantization noise σj(x)2, as shown in Appendix H.3. Thus, in the dimension with meaningful information, σj(x)
2 is much smaller than the prior variance 1, and the approximation in Eq.6 is reasonable. Let p(µ(x)) be ∏n j=1 p(µj(x)). p(µ(x)) = p(x) |det(∂x/∂z)| holds where det(∂x/∂z) is a Jacobian determinant at z = µ(x). Let CDKL be a constant n 2 log 2πe. Then, DKL(·) is reformulated by
DKL(·) ' − log ( p(µ(x)) n∏ j=1 σj(x) ) − CDKL ' − log ( p(x) ∣∣∣∣det(∂x∂z )∣∣∣∣ n∏ j=1 σj(x) ) − CDKL . (7)
Final objective form: From Eqs. 5 and 7, the objective L′x to minimise is derived as:
L′x = n∑ j=1 σj(x) 2 txzjGxxzj − β log ( p(x) ∣∣∣∣det(∂x∂z )∣∣∣∣ n∏ j=1 σj(x) ) − CDKL . (8)

heading: 3.3 PROOF OF THE HYPOTHESIS
text: Mapping VAE to implicit isometric embedding: The minimum condition of L′x at x is examined. Let x̃zj be the j-th column vector of a cofactor matrix for Jacobian matrix ∂x/∂z. Note that d log |det(∂x/∂z)|/dxzj = x̃zj/det(∂x/∂z) holds as is also used in Kato et al. (2020). Using this equation, the derivative of L′x by xzj is described by
dL′x dxzj = 2σj(x) 2Gxxzj −
β
det (∂x/∂z) x̃zj . (9)
Note that txzk · x̃zj = det(∂x/∂z) δjk holds by the cofactor’s property. Here, · denotes the dot product, and δjk denotes the Kronecker delta. By setting Eq. 9 to zero and multiplying txzk from the left, the condition to minimize L′x is derived by the next orthogonal form of xzj :
(2σj(x) 2/β) txzkGxxzj = δjk. (10)
Here, the diagonal posterior covariance is the key for orthogonality. Next, implicit latent variable y and its j-th dimensional component yj are introduced. Set yj to zero at zj = 0. The derivative between yj and zj at µj(x) is defined by
dyj dzj ∣∣∣ zj=µj(x) = √ β 2 σj(x) −1. (11)
xyj denotes ∂x/∂yj . By applying xzj = dyj/dzj xyj to Eq. 10, xyj shows the isometric property (Han & Hong, 2006; Kato et al., 2020) in the inner product space with a metric tensorGx as follows:
txyjGxxyk = δjk. (12)
Minimum entropy of implicit isometric representation: Let L′min x be a minimum of L′x at x. Dminx and Rminx denote a coding loss and KL divergence in L′minx, respectively. By applying Eqs. 10-11 and p(zj) = (dyj/dzj) p(yj) to Eqs. 5 and 7, the following equations are derived:
L′min x = Dmin x + βRmin x, where Dmin x = nβ
2 , Rmin x = − log p(y)−
n log(βπe)
2 . (13) Here, Dmin x is derived as ∑n j=1(β/2)
txyjGxxyj = nβ/2, implying each dimensional posterior variance of the implicit isometric variable is a constant β/2. In addition, exp(−L′min x/β) = p(y) exp(Const.) ∝ p(y) ' p(x) will hold in the inner product space ofGx from the isometricity. By averaging L′minx over x ∼ p(x) and approximating this average by the integration over y ∼ p(y), the global minimum L′G is derived as:
L′G = DG + βRG, where DG = nβ
2 , RG = min p(y)
( − ∫ p(y) log p(y)dy ) − n log(βπe)
2 . (14)
The term− ∫ p(y) log p(y)dy in RG is the entropy of y. Thus, the optimal implicit isometric space is derived such that the entropy of data representation is minimum in the inner product space ofGx.
When the data manifold has a disentangled property in the given metric, each yj will capture a disentangled feature with minimum entropy, as shown in Kato et al. (2020). This is analogous to PCA for Gaussian data, which gives the disentangled representation with minimum entropy in SSE. Considering the similarity to the PCA eigenvalues, the variance of yj will indicate the importance of each dimension. In the dimensions where the variance of yj is less than β/2, σj(x) = 1, µj(x) = 0, and DKLj(x) = 0 will hold. In addition, σj(x)2 txzjGxxzj will be close to 0 because this needs not to be balanced with DKLj(x). This is similar to the case in the RD theory in Eq. 4 where σzj2 is less than d, meaning no information. As a result, Eqs. 10-14 will not hold here. Thus, latent variables with variances from the largest to the n-th withDKLj(x) > 0 are sufficient for the representation and the dimensions with DKLj(x) = 0 can be ignored, allowing the reduction of the dimension n for z.
Some approximations may be slightly violated, however, our analysis still helps to understand VAE.

heading: 3.4 DISCUSSION AND RELATIONSHIP WITH PRIOR THEORETICAL STUDIES
text: First, we show β-VAE optimum as in Eq. 14 can be interpreted as the rate-distortion optimum (Eq. 4) in RD theory when the uniform distortion d in Eq. 4 is set to β/2 in the metric defined space. H(X) = − ∫ p(x) log p(x) dx denotes a differential entropy for a set x ∈ X;x ∼ p(x). For the 1-dimensional Gaussian data x ∼ N (x, 0, σ2), H(X) = 12 log(2πeσ 2) holds. Thus, Ropt in Eq. 4
is derived as a difference of the differential entropy between transformed data z ∼ ∏ j N (zj ; 0, σzj ) and uniform distortion D ∼ N (D; 0, dIm). RG is also derived as a difference of the differential entropy between transformed data y ∼ p(y) and uniform distortion D ∼ N (D; 0, (β/2)Im). Furthermore, DG in Eq. 14 can be interpreted as Dopt in Eq. 4 by setting d = β/2. As a result, the VAE optimal corresponds to the rate-distortion optimal of transform coding in RD theory, and β/2 is regarded as a variance of the constant distortion equally added to each dimensional component. Because of the isometricity, the power of distortion (i.e., posterior variance) in the implicit isometric space is the same as that in the metric defined input space. Thus the conditional distribution after optimization in the metric defined space is derived as pθ(x|z) = pθ(x|x̂) ' N (x; x̂, (β/2)I). This is consistent with the fact that the quality of the reconstructed data becomes worse in larger β.
Next, we estimate the reconstruction loss Eqφ(z|x)[log pθ(x|z)] and KL divergence DKL(·) in βVAE and also correct the analysis in Alemi et al. (2018). LetH = −Ep(x)[log p(x)] be a differential entropy of input data. When β = 1, Alemi et al. (2018) suggest ”the ELBO objective alone (and the marginal likelihood) cannot distinguish between models that make no use of the latent variable (autodecoders) versus models that make large use of the latent variable and learn useful representations for reconstruction (autoencoders),” because the reconstruction loss and KL divergence can be arbitrary value on the line−Eqφ(z|x)[log pθ(x|z)]+DKL(·) = H . Correctly, the reconstruction loss and KL divergence after optimization are deterministically estimated at any β (including β = 1) as:
Eqφ(z|x)[log pθ(x|z)] ' −(n/2) log(βπe), DKL(·) ' − log p(y)− (n/2) log(βπe). (15) The proof is explained in Appendix A.1. Thus ELBO can be estimated as:
ELBO = Ep(x)[Ez∼qφ(z|x)[log pθ(x|z)]−DKL(·)] ' Ep(x)[log p(y)] ' Ep(x)[log p(x)]. (16) As a result, when the objective of β-VAE is optimised, ELBO (Eq. 1) in the original form (Kingma & Welling, 2014) is approximately equal to the log-likelihood of x, regardless β = 1 or not.
Finally, the predetermined conditional distribution pRp(x|x̂) and the true conditional distribution after optimization pRθ(x|x̂) are examined using β in the input Euclidean space of x. Assume pRp(x|x̂) = N (x; x̂, σ2I). In this case, the metric D(x, x̂) is derived as − log pRp(x|x̂) = (1/2σ2)|x− x̂|22 + Const. From Eq. 13, the following equations are derived:
Eqφ(x̂|x)[D(x, x̂)] = Eqφ(x̂|x) [ 1 2σ2 |x− x̂|22 ] = Eqφ(x̂|x) [ 1 2σ2 ∑ i (xi − x̂i)2 ] ' nβ/2, (17)
Eqφ(x̂|x) [ (xi − x̂i)2 ] ' βσ2. (18)
Because the variance of each dimension is estimated as βσ2, the true conditional distribution after optimization is approximated as pRθ(x|x̂) ' N (x; x̂, βσ2I). If β = 1, i.e., the original VAE, pRp(x|x̂) and pRθ(x|x̂) are equivalent as expected. If β 6= 1, however, pRp(x|x̂) and pRθ(x|x̂) are different. Actually, what β-VAE does is only to scale the variance of the pre-determined conditional distribution in the original VAE by a factor of β, because β-VAE objective can be rewritten as:
Eqφ(·)[logN (x; x̂, σ 2I)]− βDKL(·) = β ( Eqφ(·)[logN (x; x̂, βσ 2I)]−DKL(·) ) + const. (19)
More detailed discussions about prior works (Higgins et al. (2017); Alemi et al. (2018); Dai et al. (2018); Dai & Wipf (2019); Tishby et al. (1999); Goyal (2001)) are explained in Appendix A.

heading: 3.5 QUANTITATIVE PROPERTIES TO VALIDATE THE THEORY
text: This section shows three quantitative properties in VAE with a priorN (z; 0, In), to validate the theory in section 3.3. The second and third properties also provide practical data analysis approaches. The derivation of equations in the second and third properties are explained in appendix C.
Norm of xyj equal to 1: Let e(j) be a vector (0, · · · , j-th
1 , · · · , 0) where the j-th dimension is 1, and others are 0. Let D′j(z) be D(Decθ(z),Decθ(z + e
(j)))/ 2, where denotes a minute value for the numerical differential. From Eq. 10, the squared norm of xyj can be numerically evaluated as the first term of Eq. 20. This value will be equal to 1 at any x and dimension j except DKLj(x) = 0.
2 β σj(x) 2D′j(z) ' 2 β
( σj(x) 2 txzjGxxzj ) ' txyjGxxyj = 1. (20)
If observed, the existence of an implicit isometric embedding can be shown because of unit norm and orthogonality (Rolı́nek et al., 2019). Eq. 20 also show σj(x)2 txzjGxxzj ' β 2 , implying that a noise σj(x) added to each dimension of latent variable causes an equal noise β/2 in the input space.
PCA-like feature: When the data manifold has a disentangled property in the given metric, the variance of the j-th implicit latent component yj can be roughly estimated as∫
yj 2p(yj)dyj '
β 2 E x∼p(x) [σj(x) −2]. (21)
The averageE[σj(x)−2] on the right allows evaluating the quantitative importance of each dimension in practice, like the eigenvalue of PCA. Note that a dimension whose average is close to 1 implies DKLj(x) = 0. Such a dimension has no information and is an exceptions of the property in Eq. 20.
Estimation of the data probability distribution: First, assume the case m = n. Since the y space is isometric to the inner product space of Gx, the PDFs in both spaces are the same. The Jacobian determinant between the input space and inner product space, giving the the ratio of PDFs, is derived as |Gx| 1 2 . We set p(µ(x)) to the prior. Thus, the data probability in the input space can be estimated by |Gx| 1 2 and either the prior/posterior or Lx after training, as the following last two equations:
p(x) ' |Gx| 1 2 p(y) ∝ |Gx| 1 2 p(µ(x)) m∏ j=1 σj(x) ∝ |Gx| 1 2 exp ( − 1 β Lx ) . (22)
In the case m > n, the derivation of the PDF ratio between the input space and the inner product space is generally intractable, except for Gx = axIm, where ax is an x-dependent scalar factor. In this case, the PDF ratio is given by axn/2. Thus, p(x) can be estimated as follows:
p(x) ∝ ax n 2 p(µ(x)) n∏ j=1 σj(x) ∝ ax n 2 exp ( − 1 β Lx ) . (23)
Equations 22 and 23 enable a probability-based quantitative data analysis/sampling in practice.

heading: 4 EXPERIMENT
text: We show the experiments of the quantitative properties presented in Section 3.5. First, the results of the toy dataset are presented. Then, the results of CelebA are shown as a real data example.

heading: 4.1 EVALUATION OF QUANTITATIVE PROPERTIES IN THE TOY DATASET
text: The toy dataset is generated as follows. First, three dimensional variables s1, s2, and s3 are sampled in accordance with the three different shapes of distributions p(s1), p(s2), and p(s3), as shown in Fig. 2. The variances of s1, s2, and s3 are 1/6, 2/3, and 8/3, respectively, such that the ratio of the variances is 1:4:16. Second, three 16-dimensional uncorrelated vectors v1, v2, and v3 with L2 norm 1, are provided. Finally, 50, 000 toy data with 16 dimensions are generated by x = ∑3 i=1 sivi. The data generation probability p(x) is also set to p(s1)p(s2)p(s3). If our hypothesis is correct, p(yj) will be close to p(sj). Then, σj(x) ∝ dzj/dyj = p(yj)/p(zj) will also vary a lot with these varieties of PDFs. Because the properties presented in Section 3.5 are calculated from σj(x), our theory can be easily validated by evaluating those properties.
Then, the VAE model is trained using Eq. 1. We use two kinds of the reconstruction loss D(·, ·) to analyze the effect of the loss metrics. The first is the square error loss equivalent to sum square error (SSE). The second is the downward-convex loss which we design as Eq. 24, such that the shape becomes similar to the BCE loss as in Appendix H.2:
D(x, x̂) = ax‖x− x̂‖22, where ax = (2/3 + 2 ‖x‖22/21) andGx = axIm. (24)
Here, ax is chosen such that the mean of ax for the toy dataset is 1.0 since the variance of x is 1/6+2/3+8/3=7/2. The details of the networks and training conditions are written in Appendix D.1.
Then the network is trained with two types of reconstruction losses. The ratio of transform loss to coding loss for the square error loss is 0.023, and that for the downward-convex loss is 0.024. As expected in section 3.2, the transform losses are negligibly small. Tables 1 and 2 show the measurements of 2βσj(x) 2D′j(z) (shown as 2 βσj 2D′j), D ′ j(z), and σj(x)
−2 described in Section 3.5. In these tables, z1, z2, and z3 show acquired latent variables. ”Av.” and ”SD” are the average and standard deviation, respectively. To begin with, the norm of the implicit orthonormal basis
Table 1: Property measurements of the toy dataset trained with the square error loss.
Table 2: Property measurements of the toy dataset trained with the downward-convex loss.
a 3/2 x p(µ(x))
∏
j σj(x), and (d) a 3/2 x exp(−Lx/β).
is discussed. In both tables, the values of 2βσ(x)j 2D′j(z) are close to 1.0 in each dimension as described in Eq. 23. By contrast, the average of D′j(z), which corresponds to txzjGxxzk , is different in each dimension. Therefore, the derivative of x with zj , the original latent variable of VAE, is not normalized.
Next, the PCA-like feature is examined. The average of σj(x)−2 in Eq.21 and its ratio are shown in Tables 1 and 2. Although the average of σj(x)−2 is a rough estimation of variance, the ratio is close to 1:4:16, i.e., the variance ratio of generation parameters s1, s2, and s3. When comparing both losses, the ratio of s2 and s3 for the downward-convex loss is somewhat smaller than that for the square error. This is explained as follows. In the downward-convex loss, |xyj |2 tends to be 1/ax from Eq. 12, i.e. txyj (axIm)xyk = δjk. Therefore, the region in the inner product space with a larger norm is shrunk, and the estimated variances corresponding to s2 and s3 become smaller.
Figure 3 shows the scattering plots of the data generation probability p(x) and estimated probabilities for the downward-convex loss. The plots for the square error loss are shown in Appendix E. Figure 3a shows the plots of p(x) and the prior probabilities p(µ(x)). This graph implies that it is difficult to estimate p(x) only from the prior. The correlation coefficient shown as ”R” (0.434) is also low. Figure 3b shows the plots of p(x) and exp(−Lx/β), i.e., the lower bound of likelihood. The correlation coefficient (0.771) becomes better, but is still not high. Next, Figures 3c and 3d show the plots of a3/2x p(µ(x)) ∏ j σj(x) and a 3/2 x exp(−Lx/β) in Eq. 23. These graphs, showing a high correlation coefficients around 0.91, support that the objective Lx in Eq. 3 is optimized in the inner product space ofGx. In the case of the square error loss, the plots with exp(−Lx/β) also shows a high correlation coefficient 0.904 because ax is 1, allowing the probability estimation from Lx in Eq. 3. The ablation study with different PDF, losses, and β is shown in Appendix E.

heading: 4.2 EVALUATIONS IN CELEBA DATASET
text: This section evaluates the first and second quantitative properties of VAE trained with the CelebA dataset 1 (Liu et al., 2015) as an example of real data. This dataset is composed of 202,599 celebrity images. In use, the images are center-cropped to form 64× 64 sized images.
1(http://mmlab.ie.cuhk.edu.hk/projects/CelebA.html)
Figure 4: Graph of σj(x)−2 average and 2βσj(x) 2D′j(z) in VAE for CelebA dataset.
Figure 5: Graph of σj(x)−2 average and 2βσj(x) 2D′j(z) in VAE for CelebA dataset with explicit decomposed loss.
We use SSIM, which is popular in image compression, as a reconstruction loss. The details of networks and training conditions are written in Appendix D.2.
Figure 4 shows the averages of σj(x)−2 in Eq.21 as the estimated variances, as well as the average and the standard deviation of 2βσj(x)
2D′j(z) in Eq.20 as the estimated square norm of implicit transform. The latent variables zi are numbered in descending order by the estimated variance. In the dimensions greater than the 27th, the averages of σj(x)−2 are close to 1 and that of 2βσj(x) 2D′j(z) is close to 0, implying DKL(·) = 0. Between the 1st and 26th dimensions, the mean and standard deviation of 2βσj(x)
2D′j(z) averages are 1.83 and 0.13, respectively. These values seem almost constant with a small standard deviation; however, the mean is somewhat larger than the expected value 1. This result implies that the implicit transform can be considered as almost orthonormal by dividing √ 1.83 ' 1.35. Thus, the average of σj(x)−2 still can determine the quantitative importance of each latent variable. This also mean that the added noise to each yj is around 1.83(β/2). We also train VAE by the decomposed loss explicitly, where Lx is set to D(x, x̆) + D(x̆, x̂) + βDKL(·). Figure 5 shows the result. Here, the mean and standard deviation of 2βσj(x)
2D′j(z) averages are 0.92 and 0.04, respectively, which suggests almost a unit norm. As a result, the explicit use of decomposed loss matches the theory better, allowing better analysis. The slight violation of the norm in the conventional form needs a more exact analysis as a future study.
Figure 6 shows decoder outputs where the selected latent variables are traversed from −2 to 2 while setting the rest to 0. The average of σj(x)−2 is also shown there. The components are grouped by the average of σj(x)−2, such that z1, z2, z3 to the large, z16, z17 to the medium, and z32 to the small, accordingly. In the large group, significant changes of background brightness, the direction of the face, and hair color are observed. In the medium group, we can see minor changes such as facial expressions. However, in the small group, there are almost no changes. This result strongly supports that the average of σj(x)−2 shows the importance of each latent variable. The traversed outputs for all the component and results with another conditions are shown in Appendix F.

heading: 5 CONCLUSION
text: This paper provides a quantitative understanding of VAE by non-linear mapping to an isometric embedding. According to the Rate-distortion theory, the optimal transform coding is achieved by using PCA/KLT orthonormal transform, where the transform space is isometric to the input. From this analogy, we show theoretically and experimentally that VAE can be mapped to an implicit isometric embedding with a scale factor derived from the posterior parameter. Based on this property, we also clarify that VAE can provide a practical quantitative analysis of input data such as the probability estimation in the input space and the PCA-like quantitative multivariate analysis. We believe the quantitative properties thoroughly uncovered in this paper will be a milestone to further advance the information theory-based generative models such as VAE in the right direction.

heading: A DETAILED RELATION TO PRIOR WORKS
text: Firstly, we clarify the the difference between ELBO in Eq. 1 and the objective Lx in Eq. 3 in the right direction. Then, we discuss the relation to the prior works. We also point out the incorrectness of some works.

heading: A.1 DERIVATION OF ELBO WITH CLEAR AND QUANTITATIVE FORM
text: We derive the reconstruction loss and KL divergence terms in ELBO (without β) at x in Eq. 1 when the objective of β-VAE L′x in Eq. 13 is optimised. The reconstruction loss can be rewritten as
Ez∼qφ(z|x)[log pθ(x|z)] = ∫ qφ(z|x) log pθ(x|z)dz = ∫ qφ(y|x) log pθ(x|y)dy. (25)
Let µy(x) be a implicit isometric variable corresponding to µ(x). Because the posterior variance in each isometric latent variable is a constant β/2, qφ(y|x) ' N (y;µy(x), (β/2)In) will hold. If β/2 is small, p(x̂) ' p(x) will hold. Then, the next equation will hold also using isometricity;
pθ(x|z) = pθ(x|y) = pθ(x|x̂) = p(x̂|x)p(x)/p(x̂) ' p(x̂|x) ' qφ(y|x). (26)
Thus the reconstruction loss is estimated as: Ez∼qφ(z|x)[log pθ(x|z)] ∼ ∫ N (y;µy(x), (β/2)In) logN (y;µy(x), (β/2)In) dy
= −(n/2) log(βπe). (27)
From Eq. 13, KL divergence is derived as:
DKL(·) = Rminx = − log p(y)− (n/2) log(βπe). (28)
By summing both terms, ELBO at x can be estimated as
ELBO = Ex∼p(x)[Ez∼qφ(z|x)[log pθ(x|z)]−DKL(·)] ' Ex∼p(x)[log p(y)] ' Ex∼p(x)[log p(x)]. (29)
As a result, ELBO (Eq. 1) in the original form (Kingma & Welling, 2014) is close to the loglikelihood of x, regardless β = 1 or not, when the objective of β-VAE (Higgins et al., 2017) is optimised.
Some of the prior VAE works do not explicitly distinguish between the reconstruction loss Ex̂∼p(x̂|x)[p(x|x̂)] in ELBO and the distortion D(x, x̂) in the objective Lx by mistake, which leads to some incorrect discussion.
In addition, there have also been incorrect discussions in some prior works. In ELBO derivation, they useEx̂∼p(x̂|x)[p(x|x̂)] as a reconstruction loss, without discussing what kinds of properties the distortion probability p(x̂|x) should be. In training VAE with a real dataset, by contrast, they use a predetermined distortion metric D(x, x̂) like BCE and SSE as a reconstruction loss instead of a log-likelihood of the distortion probability, without discussing what distortion probability should be after optimization.
Correctly, the distortion probability p(x̂|x) after training is determined by β and the metric as pθ(x̂|x) ' N (x̂;x, (β/2)Im) in the metric defined space. Then, by applying pθ(x̂|x) ' pθ(x̂|x) to Eq. 1, the value of ELBO will become log p(y) ' log p(x) regardless β = 1 or not.
IfD(x,x+δx) = tδxGxδx+O(||δx||3) is not SSE, by introducing a variable x́ = Lx−1x where Lx satisfies tLx Lx = Gx, the metric D(·, ·) can be replaced by SSE in the Euclidean space of x́.

heading: A.2 RELATION TO TISHBY ET AL. (1999)
text: The theory described in Tishby et al. (1999) is consistent with our analysis. Tishby et al. (1999) clarified the behaviour of the compressed representation when the rate-distortion trade-off is optimized.
x ∈ X denotes the signal space with a fixed probability p(x) and x̂ ∈ X̂ denotes its compressed representation. Let D(x, x̂) be a loss metric. Then the rate-distortion trade-off can be described as:
L = I(X; X̂) + β′ E p(x,x̂) [D(x, x̂)]. (30)
By solving this condition, they derive the following equation:
p(x̂|x) ∝ exp(−β′D(x, x̂)). (31)
As shown in our discussion above, p(x̂|x) ' N (x̂;x, (β/2)Im) will hold in the metric defined space from our VAE analysis. This result is equivalent to Eq. 31 in their work if D(x, x̂) is SSE and β′ is set to β−1, as follows:
p(x̂|x) ∝ exp(−β′D(x, x̂)) = exp ( −||x− x̂|| 2 2
2(β/2)
) ∝ N (x̂;x, (β/2)Im). (32)
If D(x, x̂) is not SSE, the use of the space transformation explained in appendix A.1 will lead to the same result.
A.3 RELATION TO β-VAE (HIGGINS ET AL., 2017)
In the β-VAE work by Higgins et al. (2017), it is presumed that the objective Lx was not mistakenly distinguished from ELBO. In their work, ELBO equation is modified as:
Ep(x)[ Ex̂∼pφ(x̂|x)[qθ(x|x̂)]− βDKL(·) ]. (33)
However, they use the predetermined probabilities of pθ(x̂|x) such as the Bernoulli and Gaussian distributions in training (described in table 1 in Higgins et al. (2017)). As shown in our appendix H.2, the log-likelihoods of the Bernoulli and Gaussian distributions can be regarded as BCE and SSE metrics, respectively. As a result, the actual objective for training in Higgins et al. (2017) is not Eq. 33, but the objective Lx in Eq. 3 using BCE and SSE metrics with varying β. Thus ELBO as Eq. 1 form will become log p(x) in the BCE / SSE metric defined space regardless β = 1 or not, as shown in appendix A.1.
Actually, the equation 33 dose not show the log-likelihood of x. When DKL(·) ' − log p(x) − (n/2) log(βπe) and Ex̂∼p(x̂|x)[p(x|x̂)] ' −(n/2) log(βπe) are applied, the value of Eq. 33 is derived as β log p(x) + (β − 1)(n/2) log(βπe), which is different from the log-likelihood of x if β 6= 1. Correctly, what β-VAE really does is only to scale the variance of the pre-determined conditional distribution in the original VAE by a factor of β. In the case the pre-determined conditional distribution is Gaussian N (x; x̂, σ2I), the objective of β can be can be rewritten as a linearly scaled original VAE objective with a Gaussian N (x; x̂, βσ2I):
Eqφ(·)[logN (x; x̂, σ 2I)]− βDKL(·) = Eqφ(·)
[ −1
2 log 2πσ2 − |x− x̂| 2 2 2σ2
] − βDKL(·)
= β ( Eqφ(·) [ −1 2 log 2πβσ2 − |x− x̂| 2 2 2βσ2 ] −DKL(·) ) + β
2 log 2πβσ2 − 1 2 log 2πσ2
= β ( Eqφ(·)[logN (x; x̂, βσ 2I)]−DKL(·) ) + const.
(34)

heading: A.4 RELATION TO ALEMI ET AL. (2018)
text: Alemi et al. (2018) discuss the rate-distortion trade-off by the theoretical entropy analysis. Their work is also presumed that the objective Lx was not mistakenly distinguished from ELBO, which leads to the incorrect discussion. In their work, the differential entropy for the inputH , distortionD, and rate R are derived carefully. They suggest that VAE with β = 1 is sensitive (unstable) because D and R can be arbitrary value on the line R = H − βD = H −D. Furthermore, they also suggest
that R ≥ H, D = 0 at β → 0 and R = 0, D ≥ H at β → ∞ will hold as shown the figure 1 of their work.
In this appendix, we will show that β determines the value of R and D specifically. We also show that R ' H −D will hold regardless β = 1 or not. In their work, these values of H , D, and D are mathematically defined as:
H ≡ − ∫ dx p∗(x) log p∗(x), (35)
D ≡ − ∫ dx p∗(x) ∫ dz e(z|x) log d(x|z), (36)
R ≡ ∫ dx p∗(x) ∫ dz e(z|x) log e(z|x)
m(z) . (37)
Here, p∗(x) is a true PDF of x, e(z|x) is a stochastic encoder, e(z|x) is a decoder, and m(z) is a marginal probability of z.
Our work allows a rough estimation of Eqs. 35-37 with β by introducing the implicit isometric variable y as explained in our work.
Using isometric variable y and the relation dz e(z|x) = dy e(y|x), Eq. 36 can be rewritten as: D = − ∫ dx p∗(x) ∫ dy e(y|x) log d(x|y). (38)
Let µy be the implicit isometric latent variable corresponding to the mean of encoder output µ(x). As discussed in section 3.3, e(y|x) = N (y;µy, (β/2)In) will hold. Because of isometricity, the value of d(x|y) will be also close to e(y|x) = N (y;µy, (β/2)In). Though d(x|z) must depend on e(z|x), this important point has not been discussed well in this work. By using the implicit isometric variable, we can connect both theoretically. Thus, D can be estimated as:
D ' ∫ dx p∗(x) ∫ dy N (y;µy, (β/2)In) logN (y;µy, (β/2)In)
' ∫ dx p∗(x) (n
2 log(βπe) ) = n
2 log(βπe) (39)
Second, R is examined. m(y) is a marginal probability of y. Using the relation dz e(z|x) = dy e(y|x) and e(z|x)/m(z) = (e(y|x)(dy/dz))/(m(y)(dy/dz)) = e(y|x)/m(y), Eq. 37 can be rewritten as:
R ' ∫ dx p∗(x) ∫ dy e(y|x) log e(y|x)
m(y) . (40)
Because of isometricity, e(y|x) ' p(x̂|x) ' N (x̂;x, (β/2)Im) will approximately hold where x̂ denotes a decoder output. Thus m(y) can be approximated by:
m(y) ' ∫ dx p∗(x)e(y|x) ' ∫ dx p∗(x) N (x̂;x, (β/2)Im) (41)
Here, if β/2, i.e., added noise, is small enough compared to the variance of x, a normal distribution function term in this equation will act like a delta function. Thus m(y) can be approximated as:
m(y) ' ∫ dx́ p∗(x́) δ(x́− x) ' p∗(x). (42)
In the similar way, the following approximation will also hold.∫ dy e(y|x) logm(y) ' ∫ dy e(y|x) log p∗(x) ' ∫ dx́ δ(x́− x) log p∗(x́) ' log p∗(x) (43)
By using these approximation and applying Eqs. 38-39, R in Eq. 37 can be approximated as: R ' ∫ dx p∗(x) ∫ dy e(y|x) log e(y|x)
p∗(x) ' − ∫ dx p∗(x) log p∗(x)− ( − ∫ dx p∗(x) ∫ dy e(y|x) log e(y|x) ) ' H − n
2 log(βπe)
' H −D (44) As discussed above, R and D can be specifically derived from β. In addition, Shannon lower bound discussed in Alemi et al. (2018) can be roughly verified in the optimized VAE with clearer notations using β.
From the discussion above, we presume Alemi et al. (2018) might wrongly treat D in their work. They suggest that VAE with β = 1 is sensitive (unstable) becauseD andR can be arbitrary value on the line R = H − βD = H −D; however, our work as well as Tishby et al. (1999) (appendix A.2) and Dai & Wipf (2019)(appendix A.5) show that the differential entropy of the distortion and rate, i.e., D and R, are specifically determined by β after optimization, and R = H − D will hold for any β regardless β = 1 or not. Alemi et al. (2018) also suggest D should satisfy D ≥ 0 because D is a distortion; however, we suggest D should be treated as a differential entropy and can be less than 0 because x is once handled as a continuous signal with a stochastic process in Eqs. 35-37. Here, D ' (n/2) log(βπe) can be −∞ if β → 0, as also shown in Dai & Wipf (2019). Thus, upper bound of R at β → 0 is not H , but R = H − (−∞) =∞, as shown in RD theory for a continuous signal. Huang et al. (2020) show this property experimentally in their figures 4-8 such that R seems to diverge if MSE is close to 0.

heading: A.5 RELATION TO DAI ET AL. (2018) AND DAI & WIPF (2019)
text: Our work is consistent with Dai et al. (2018) and Dai & Wipf (2019).
Dai et al. (2018) analyses VAE by assuming a linear model. As a result, the estimated posterior is constant. If the distribution of the manifold is the Gaussian, our work and Dai et al. (2018) give a similar result with constant posterior variances. For non-Gaussian data, however, the quantitative analysis such as probability estimation is intractable using their linear model. Our work reveals that the posterior variance gives a scaling factor between z in VAE and y in the isometric space when VAE is ideally trained with rich parameters. This is validated by Figures 3c and 3d, where the estimation of the posterior variance at each data point is a key.
Next, the relation to Dai & Wipf (2019) is discussed. They analyse a behavior of VAE when ideally trained. For example, the theorem 5 in their work shows that D → (d/2) log γ + O(1) and R → −(γ̂/2) log γ+O(1) hold if γ → +0, where γ, d, and γ̂ denote a variance of d(x|z), data dimension, and latent dimension, respectively. By setting γ = β/2 and d = γ̂ = n, this theorem is consistent with R and D derived in Eq. 39 and Eq. 44.

heading: A.6 RELATION TO TRANSFORM CODING (GOYAL, 2001)
text: We show the optimum condition of VAE shown in Eq. 14 can be mapped to the optimum condition of transform coding (Goyal, 2001) as shown in Eq. 4. First, the derivation of Eq. 4 is explained by solving the optimal distortion assignment to each dimension. In the transform coding for m dimensional the Gaussian data, an input data x is transformed to z using an orthonormal transform such as KLT/DCT. Then each dimensional component zj is encoded with allowing distortion dj . Let D be a target distortion satisfying D = ∑m j=1 dj . σ 2 zj denotes a variance of each dimensional
component zj for the input dataset. Then, a rate R can be derived as ∑m j=1 1 2 log(σ 2 zj/dj). By introducing a Lagrange parameter λ and minimizing a rate-distortion optimization costL = D+λR, the optimum condition is derived as:
λopt = 2D/m, dj = D/m = λopt/2. (45) This result is consistent with Eq. 14 by setting β = λopt = 2D/m. This implies that L′G in Eq. 14 is a rate-distortion optimization (RDO) cost of transform coding when x is deterministically transformed to y in the implicit isometric space and stochastically encoded with a distortion β/2.
B ESTIMATION OF THE CODING LOSS AND TRANSFORM LOSS IN 1-DIMENSIONAL LINEAR VAE
This appendix estimates the coding loss and transform loss in 1-dimensional linear β-VAE for the Gaussian data, and also shows that the result is consistent with the Wiener filter. Let x be a one dimensional data with the normal distribution:
x ∈ R, x ∼ N (x; 0, σx2) (46)
Let z be a one dimensional latent variable. Following two linear encoder and decoder are provided with constant parameters a, b, and σz to optimize:
z = ax+ σz where ∼ N ( ; 0, 1), x̂ = bz. (47)
First, KL divergence at x, DKLx is derived. Due to the above relationship, we have
p(z) = N (z; 0, (aσx)2). (48)
Using Eq. 6, KL-divergence at x can be evaluated as:
DKLx = − log(σzp(z))− 1
2 log 2πe = − log σx +
a2x2
2 − 1 2 . (49)
Second, the reconstruction loss at x Dx is evaluated as:
Dx = E ∼N ( ;0,1)
[(x− (b(ax+ σz )))2] = ((ab− 1)x)2 + b2σz2. (50)
Then, the loss objective Lx = Dx + βDKLx is averaged over x ∼ N (x; 0, σx2), and the objective L to minimize is derived as:
L = E x∼N (x;0,σx2)
[Lx] = (ab− 1)2σx2 + b2σz2 + β ( − log σz + a2σx 2
2 − 1 2
) . (51)
Here, (ab − 1)2σx2 and b2σz2 in the last equation are corresponding to the transform loss DT and coding loss DC, respectively.
By solving dL/da = 0, dL/db = 0, and dL/dσz = 0, a , b, and σz are derived as follows:
a = 1/σx,
b = σx
( 1 + √ 1− 2β/σx2 ) 2 ,
σz = 2 √ β/2
σx
( 1 + √ 1− 2β/σx2 ) . (52) From Eq. 52, DT and DC are derived as:
DT =
(√ 1− 2β/σx2 − 1
2
)2 σx 2,
DC = β/2. (53)
As shown in section 3.3, the added noise, β/2, should be reasonably smaller than the data variance σx 2. If σx2 β, b and σz in Eq. 52 can be approximated as:
DT ' (β/2)2 σx2 = β/2 σx2 DC. (54)
As shown in this equation, DT/DC is small in the VAE where the added noise is reasonably small, and DT can be ignored.
Next, the relation to the Wiener filter is discussed. We consider an simple 1-dimensional Gaussian process. Let x ∼ N (x; 0, σ2x) be input data. Then, x is scaled by s, and a Gaussian noise n ∼
N (n; 0, σ2n) is added. Thus, y = s x + n is observed. From the Wiener filter theory, the estimated value with minimum distortion, x̂ can be formulated as:
x̂ = sσx
2
s2σx2 + σn2 y. (55)
In this case, the estimation error is derived as:
E[(x̂− x)2] = σn 4
(s2σx2 + σn2)2 σx
2 + s2σx 4
(s2σx2 + σn2)2 σn
2 = σx
2
σx2 + (σn2/s2) (σn
2/s2). (56)
In the second equation, the first term is corresponding to the transform loss, and the second term is corresponding to the coding loss. Here the ratio of the transform loss and coding loss is derived as σn 2/(s2σx 2). By appying s = 1/σx and σn = σz to σn2/(s2σx2) and assuming σ2x β/2, this ratio can be described as: σn 2
s2σx2 = σz
2 = β/2
σ2x 4( 1 + √ 1− 2β/σx2 )2 = β/2σ2x +O (( β/2 σ2x )2) . (57)
This result is consistent with Eq. 54, implying that optimized VAE and the Wiener filter show similar behaviours.

heading: C DERIVATION OF QUANTITATIVE PROPERTIES IN SECTION 3.5

heading: C.1 DERIVATION OF THE ESTIMATED VARIANCE
text: This appendix explains the derivation of Eq. 21 in Section 3.5. Here, we assume that zj is mapped to yj such that yj is set to 0 at zj = 0. We also assume that the prior distribution is N (z; 0, In). The variance is derived by the subtraction of E[yj ]
2, the square of the mean, from E[y2j ], the square mean. Thus, the approximations of both E[yj ] and E[y2j ] are needed.
First, the approximation of the mean E[yj ] is explained. Because the cumulative distribution functions (CDFs) of yj are the same as CDF of zj , the following equations hold:∫ 0
−∞ p(yj)dyj = ∫ 0 −∞ p(zj)dzj = 0.5, ∫ ∞ 0 p(yj)dyj = ∫ ∞ 0 p(zj)dzj = 0.5. (58)
This equation means that the median of the yj distribution is 0. Because the mean and median are close in most cases, the mean E[yj ] can be approximated as 0. As a result, the variance of yj can be approximated by the square mean E[y2j ].
Second, the approximation of the square mean E[y2j ] is explained. The standard deviation of the posterior σj(x) is assumed as a function of zj , regardless of x. This function is denoted as σj(zj). For zj ≥ 0, yj is approximated as follows, using Eq. 11 and replacing the average of 1/σj(źj) over źj = [0, zj ] by 1/σj(zj):
yj = ∫ zj 0 dyj dźj dźj = √ β 2 ∫ zi 0 1 σj(źj) dźi ' √ β 2 1 σj(zj) ∫ zj 0 dźj = √ β 2 zj σj(zj) . (59)
The same approximation is applied to zi < 0. Then the square mean of yi is approximated as follows, assuming that the correlation between σ(zj)
−2 and zj2 is low:∫ yj 2p(yj)dyj ' β
2
∫ ( zj
σj(zj)
)2 p(zj)dzj ' β
2
∫ σj(zj) −2 p(zj)dzj ∫ zj 2p(zj)dzj . (60)
Finally, the square mean of yi is approximated as the following equation, using ∫ zj 2p(zj)dzj = 1 and replacing σj(zj) 2 by σj(x)2, i.e., the posterior variance derived from the input data:∫
yj 2p(yj)dyj '
β
2
∫ σj(zj) −2 p(zj)dzj ' β
2 E zj∼p(zj) [σj(zj)
−2 ] ' β
2 E x∼p(x) [σj(x)
−2]. (61)
Although some rough approximations are used in the expansion, the estimated variance in the last equation seems still reasonable, because σj(x) shows a scale factor between yj and zj while the variance of zj is always 1 for the priorN (zj ; 0, 1). Considering the variance of the prior ∫ zj
2p(zj)dzj in the expansion, this estimation method can be applied to any prior distribution.

heading: C.2 DERIVATION OF THE DATA PROBABILITY ESTIMATION
text: This appendix shows the derivation of variables in Eqs. 22 and 23. First, the derivation of Lx for the input x is described. Then, the PDF ratio between the input space and inner product space is explained for the cases m = n and m > n.

heading: Derivation of Lx for the input x :
text: As shown in in Eq. 1,Lx is denoted as−Ez∼qφ(z|x)[ · ]+βDKL( · ). We approximateEz∼qφ(z|x)[ · ] as 12 (D(x,Decθ(µx +σx)) +D(x,Decθ(µx −σx))), i.e., the average of two samples, instead of the average over z ∼ qφ(z|x). DKL( · ) can be calculated from µx and σx using Eq. 2.

heading: The PDF ratio in the case m = n:
text: The PDF ratio form = n is a Jacobian determinant between two spaces. First, (∂x∂y ) TGx( ∂x ∂y ) = Im holds from Eq. 12. |∂x/∂y|2 |Gx| = 1 also holds by calculating the determinant. Finally, |∂x/∂y| is derived as |Gx|1/2 using |∂y/∂x| = |∂x/∂y|−1. The PDF ratio in the case m > n andGx = axIm: Although the strict derivation needs the treatment of the Riemannian manifold, we provide a simple explanation in this appendix. Here, it is assumed that DKL(j)(·) > 0 holds for all j = [1, ..n]. If DKL(j)(·) = 0 for some j, n is replaced by the number of latent variables with DKL(j)(·) > 0.
For the implicit isometric space Siso(⊂ Rm), there exists a matrix Lx such that both y = Lxx and Gx =
tLxLx holds. w denotes a point in Siso, i.e., w ∈ Siso. Because Gx is assumed as axIm in Section 3.5, Lx = ax1/2Im holds. Then, the mapping function w = h(x) between Sinput and Siso is defined, such that:
∂h(x) ∂x = ∂w ∂x = Lx, and h(x(0)) = w(0) for ∃ x(0) ∈ Sinput and ∃ w(0) ∈ Siso. (62)
Let δx and δw are infinitesimal displacements around x and w = h(x), such that w + δw = h(x+ δx). Then the next equation holds from Eq. 62:
δw = Lxδx (63)
Let δx(1), δx(2), δw(1), and δw(2) be two arbitrary infinitesimal displacements around x and w = h(x), such that δw(1) = Lxδx(1) and δw(2) = Lxδx(2). Then the following equation holds, where · denotes the dot product.
tδx(1)Gxδx (2) = t(Lxδx (1))(Lxδx (2)) = δw(1) · δw(2) (64)
This equation shows the isometric mapping from the inner product space for x ∈ Sinput with the metric tensorGx to the Euclidean space for w ∈ Siso. Note that all of the column vectors in the Jacobian matrix ∂x/∂y also have a unit norm and are orthogonal to each other in the metric space for x ∈ Sinput with the metric tensor Gx. Therefore, the m×n Jacobian matrix ∂w/∂y should have a property that all of the column vectors have a unit norm and are orthogonal to each other in the Euclidean space.
Then n-dimensional space which is composed of the meaningful dimensions from the implicit isometric space is named as the implicit orthonormal space Sortho. Figure 7 shows the projection of the volume element from the implicit orthonormal space to the isometric space and input space. Let dVortho be an infinitesimal n-dimensional volume element in Sortho. This volume element is a n-dimensional rectangular solid having each edge length dyj . Let Vn(dVX) be the n-dimensional volume of a volume element dVX. Then, Vn(dVortho) = ∏n j dyj holds. Next, dVortho is projected to n dimensional infinitesimal element dViso in Siso by ∂w/∂y. Because of the orthonormality, dViso is equivalent to the rotation / reflection of dVortho, and Vn(dViso) is the same as Vn(dVortho), i.e., ∏n j dyj . Then, dViso is projected to n-dimensional element dVinput in Sinput by ∂x/∂w = L−1x = ax −1/2Im. Because each dimension is scaled equally by the scale factor ax −1/2, Vn(dVinput) = ∏n j ax −1/2dyj = ax −n/2 Vn(dVortho) holds. Here, the ratio of the volume element between Sinput and Sortho is Vn(dVinput)/Vn(dVortho) = ax−n/2. Note that the PDF ratio is derived by the reciprocal of Vn(dVinput)/Vn(dVortho). As a result, the PDF ratio is derived as axn/2.

heading: D DETAILS OF THE NETWORKS AND TRAINING CONDITIONS IN THE EXPERIMENTS
text: This appendix explains the networks and training conditions in Section 4.

heading: D.1 TOY DATA SET
text: This appendix explains the details of the networks and training conditions in the experiment of the toy data set in Section 4.1.

heading: Network configurations:
text: FC(i, o, f) denotes a FC layer with input dimension i, output dimension o, and activate function f.
The encoder network is composed of FC(16, 128, tanh)-FC(128, 64, tahh)-FC(64, 3, linear)×2 (for µ and σ). The decoder network is composed of FC(3, 64, tanh)-FC(64, 128, tahh)-FC(128, 16, linear).

heading: Training conditions:
text: The reconstruction loss D(·, ·) is derived such that the loss per input dimension is calculated and all of the losses are averaged by the input dimension m = 16. The KL divergence is derived as a summation of DKL(j)(·) as explained in Eq. 2. In our code, we use essentially the same, but a constant factor scaled loss objective from the original β-VAE form Lx = D(·, ·) + βDKL(j)(·) in Eq. 1, such as:
Lx = λ D(·, ·) +DKL(j)(·). (65)
Equation 65 is essentially equivalent to L = D(·, ·) + βDKL(j)(·), multiplying a constant λ = β−1 to the original form. The reason why we use this form is as follows. Let ELBOtrue be the true ELBO in the sense of log-likelihood, such as E[log p(x)]. As shown in Section 3.3, the minimum of the loss objective in the original β-VAE form is likely to be a −βELBOtrue + Constant. If we use Eq. 65, the minimum of the loss objective will be −ELBOtrue + Constant, which seems more natural form of ELBO. Thus, Eq. 65 allows estimating a data probability from Lx in Eqs. 22 and 23, without scaling Lx by 1/β.
Then the network is trained with λ = β−1 = 100 using 500 epochs with a batch size of 128. Here, Adam optimizer is used with the learning rate of 1e-3. We use a PC with CPU Inter(R) Xeon(R) CPU E3-1280v5@3.70GHz, 32GB memory equipped with NVIDIA GeForce GTX 1080. The simulation time for each trial is about 20 minutes, including the statistics evaluation codes.
In our experiments, λ or β−1, i.e., 100, seems somewhat large. This is caused by the use of the mean square error as a reconstruction loss. In contrast, KL divergence is the sum for the whole image, which can be thought of as a rate for the whole image. Considering the number of input dimensions, β′ = (λ/16)−1 = 16/λ = 0.16 is thought of as β in the general form of VAE.

heading: D.2 CELEBA DATA SET
text: This appendix explains the details of the networks and training conditions in the experiment of the toy data set in Section 4.2.

heading: Network configurations:
text: CNN(w, h, s, c, f) denotes a CNN layer with kernel size (w, h), stride size s, dimension c, and activate function f. GDN and IGDN 2 are activation functions designed for image compression (Ballé et al., 2016). This activation function is effective and popular in deep image compression studies.
The encoder network is composed of CNN(9, 9, 2, 64, GDN) - CNN(5, 5, 2, 64, GDN) - CNN(5, 5, 2, 64, GDN) - CNN(5, 5, 2, 64, GDN) - FC(1024, 1024, softplus) - FC(1024, 32, None)×2 (for µ and σ) in encoder.
The decoder network is composed of FC(32, 1024, softplus) - FC(1024, 1024, softplus) - CNN(5, 5, 2, 64, IGDN) - CNN(5, 5, 2, 64, IGDN) - CNN(5, 5, 2, 64, IGDN)-CNN(9, 9, 2, 3, IGDN).

heading: Training conditions:
text: In this experiment, SSIM explained in Appendix H.2 is used as a reconstruction loss. The reconstruction loss D(·, ·) is derived as follows. Let SSIM be a SSIM calculated from two input images. Then 1 − SSIM is set to D(·, ·). The KL divergence is derived as a summation of DKL(j)(·) as explained in Eq. 2.
We also use the loss form as in Equation 65 in our code. In the case of the decomposed loss, the loss function Lx is set to λ(D(x, x̆) +D(x̆, x̂)) +DKL(·) in our code. Then, the network is trained with λ = β−1 = 1, 000 using a batch size of 64 for 300,000 iterations. Here, Adam optimizer is used with the learning rate of 1e-3.
We use a PC with CPU Intel(R) Core(TM) i7-6850K CPU @ 3.60GHz, 12GB memory equipped with NVIDIA GeForce GTX 1080. The simulation time for each trial is about 180 minutes, including the statistics evaluation codes.
In our experiments, λ = β−1 = 1, 000 seems large. This is caused by the use of SSIM. As explained in Appendix H.2, SSIM is measured for a whole image, and its range is between 0 and 1. The order of 1 − SSIM is almost equivalent to the mean square error per pixel, as shown in Eq. 74. As explained in Appendix D.1, KL divergence is thought of as a rate for the whole image. Considering the number of pixels in a image, β′ = (λ/(64 × 64))−1 = 4096/λ = 4.096 is comparable to β in the general form of VAE.
2Google provides a code in the official Tensorflow library (https://github.com/tensorflow/compression)

heading: E ADDITIONAL RESULTS IN THE TOY DATASETS

heading: E.1 SCATTERING PLOTS FOR THE SQUARE ERROR LOSS IN SECTION
text: 4.1 Figure 8a shows the plots of p(x) and estimated probabilities for the square error coding loss in Section 4.1, where the scale factor ax in Eq. 23 is 1. Thus, both exp(−Lx/β) and p(µ(x)) ∏ j σj(x) show a high correlation, allowing easy estimation of the data probability in the input space. In contrast, p(µ(x)) still shows a low correlation. These results are consistent with our theory.
E.2 ABLATION STUDY USING 3 TOY DATASETS, 3 CODING LOSSES, AND 10 β PARAMETERS.
In this appendix, we explain the ablation study for the toy datasets. We introduce three toy datasets and three coding losses including those used in Section 4.1. We also change β−1 = λ from 1 to 1, 000 in training. The details of the experimental conditions are shown as follows.
Datasets: First, we call the toy dataset used in Section 4.1 the Mix dataset in order to distinguish three datasets. The second dataset is generated such that three dimensional variables s1, s2, and s3 are sampled in accordance with the distributions p(s1), p(s2), and p(s3) in Figure 9. The variances of the variables are the same as those of the Mix dataset, i.e., 1/6, 2/3, and 8/3, respectively. We call this the Ramp dataset. Because the PDF shape of this dataset is quite different from the prior N (z; 0, I3), the fitting will be the most difficult among the three. The third dataset is generated such that three dimensional variables s1, s2, and s3 are sampled in accordance with the normal distributions N (s1; 0, 1/6), N (s2; 0, 2/3), and N (s3; 0, 8/3), respectively. We call this the Norm dataset. The fitting will be the easiest, because both the prior and input have the normal distributions, and the posterior standard deviation, given by the PDF ratio at the same CDF, can be a constant.
Coding losses: Two of the three coding losses is the square error loss and the downward-convex loss described in Section 4.1. The third coding loss is a upward-convex loss which we design as Eq. 66 such that the scale factor ax becomes the reciprocal of the scale factor in Eq. 24:
D(x, x̂) = ax‖x− x̂‖22, where ax = (2/3 + 2 ‖x‖22/21)−1 andGx = axIm. (66) Figure 10 shows the scale factors ax in Eqs. 24 and 66, where s1 in x = (s1, 0, 0) moves within ±5. Parameters: As explained in Appendix D.1, λ = 1/β is used as a hyper parameter. Specifically, λ = 1, 2, 5, 10, 20, 50, 100, 200, 500, and 1, 000 are used.
Figures 11 - 19 show the property measurements for all combinations of the datasets and coding losses, with changing λ. In each Figure, the estimated norms of the implicit transform are shown in the figure (a), the ratios of the estimated variances are shown in the figure (b), and the correlation coefficients between p(x) and estimated data probabilities are shown in the figure (c), respectively.
First, the estimated norm of the implicit transform in the figures (a) is discussed. In all conditions, the norms are close to 1 as described in Eq. 20 in the λ range 50 to 1000. These results show consistency with our theoretical analysis, supporting the existence of the implicit orthonormal transform. The values in the Norm dataset are the closest to 1, and those in the Ramp dataset are the most different, which seems consistent with the difficulty of the fitting.
Second, the ratio of the estimated variances is discussed. In the figures (b), Var(zj) denotes the estimated variance, given by the average of σ−2j(x). Then, Var(z2)/Var(z1) and Var(z3)/Var(z1) are plotted. In all conditions, the ratios of Var(z2)/Var(z1) and Var(z3)/Var(z1) are close to the variance ratios of the input variables, i.e., 4 and 16, in the λ range 5 to 500. Figure 20 shows the detailed comparison of the ratio for the three datasets and three coding losses at λ = 100. In most cases, the estimated variances in the downward-convex loss are the smallest, and those in the upward-convex loss are the largest, which is more distinct for Var(z3)/Var(z1). This can be explained as follows. When using the downward-convex loss, the space region with a large norm is thought of as shrinking in the inner product space, as described in Section 4.1. This will make the variance smaller. In contrast, when using the upward-convex loss, the space region with a large norm is thought of as expanding in the inner product space, making the variance larger. Here, the dependency of the losses on the ratio changes is less in the Norm dataset. The possible reason is that data in the normal distribution concentrate around the center, having less effect on the loss scale factor in the downward-convex loss and upward-convex loss.
Third, the correlation coefficients between p(x) and the estimated data probabilities in the figures (c) are discussed. In the Mix dataset and Ramp dataset, the correlation coefficients are around 0.9 in the λ range from 20 to 200 when the estimated probabilities axn/2p(µ(x)) ∏n j=1 σj(x) and ax n/2 exp(−(1/β)Lx) in Eq. 23 are used. When using p(µ(x)) ∏n j=1 σj(x) and exp(−(1/β)Lx) in the downward-convex loss and upward-convex loss, the correlation coefficients become worse. In addition, when using the prior probability p(µ(x)), the correlation coefficients always show the worst. In the Norm dataset, the correlation coefficients are close to 1.0 in the wider range of λ when using the estimated distribution in Eq. 23. When using p(µ(x)) ∏n j=1 σj(x) and exp(−(1/β)Lx) in the downward-convex loss and upward-convex loss, the correlation coefficients also become worse. When using the prior probability p(µ(x)), however, the correlation coefficients are close to 1 in contrast to the other two datasets. This can be explained because both the input distribution and the prior distribution are the same normal distribution, allowing the posterior variances almost constant. These results also show consistency with our theoretical analysis.
Figure 21 shows the dependency of the coding loss on β for the Mix, Ramp, and Norm dataset using square the error loss. From DG in Eq. 14 and n = 3, the theoretical value of coding loss is 3β2 , as also shown in the figure. Unlike Figs. 11-19, x-axis is β = λ−1 to evaluate the linearity. As expected in section 3.3, the coding losses are close to the theoretical value where β < 0.1, i.e., λ > 10.
Figure 22 shows the dependency of the ratio of transform loss to coding loss on β for the Mix, Ramp, and Norm dataset using square the error loss. From Eq. 54, the estimated transform loss is∑3 i=1(β/2) 2/Var(si) = 63β2 32 . Thus the theoretical value is ( 63β2 32 )/( 3β 2 ) = 21β 16 , as is also shown in the figure. x-axis is also β = λ−1 like Figure 21. Considering the correlation coefficient discussed above, the useful range of β seems between 0.005-0.05 (20-200 for λ). In this range, the ratio is less than 0.1, implying the transform loss is almost negligible. As expected in section 3.2 and appendix B, the ratio is close to the theoretical value where β > 0.01, i.e., λ < 100. For β < 0.01, the conversion loss is still negligibly small, but the ratio is somewhat off the theoretical value. The reason is presumably that the transform loss is too small to fit the network.
As shown above, this ablation study strongly supports our theoretical analysis in sections 3.

heading: F.1 TRAVERSED OUTPUTS FOR ALL THE COMPONENT IN THE EXPERIMENTAL SECTION 4.2

heading: F ADDITIONAL RESULTS IN CELEBA DATASET
text: Figure 23 shows decoder outputs for all the components, where each latent variable is traversed from −2 to 2. The estimated variance of each yj , i.e., σ−2j , is also shown in these figures. The latent variables zi are numbered in descending order by the estimated variances. Figure 23a is a result using the conventional loss form, i.e., Lx = D(x, x̂) + βDKL(·). The degrees of change seem to descend in accordance with the estimated variances. In the range where j is 1 from 10, the degrees of changes are large. In the range j > 10, the degrees of changes becomes gradually smaller. Furthermore, almost no change is observed in the range j > 27. As shown in Figure 4, DKL(j)(·) is close to zero for j > 27, meaning no information. Thus, this result is clearly consistent with our theoretical analysis in Section
Figure 23b is a result using the decomposed loss form, i.e., Lx = D(x, x̆) +D(x̆, x̂) + βDKL(·). The degrees of change also seem to descend in accordance with the estimated variances. When looking at the detail, there are still minor changes even j = 32. As shown in Figure 5, KL divergences DKL(j)(·) for all the components are larger than zero. This implies all of the dimensional components have meaningful information. Therefore, we can see a minor change even j = 32. Thus, this result is also consistent with our theoretical analysis.
Another minor difference is sharpness. Although the quantitative comparison is difficult, the decoded images in Figure 23b seems somewhat sharper than those in Figure 23a. A possible reason for this minor difference is as follows. The transform loss D(x, x̆) serves to bring the decoded image of µ(x) closer to the input. In the conventional image coding, the orthonormal transform and its inverse transform are used for encoding and decoding, respectively. Therefore, the input and the decoded output are equivalent when not using quantization. If not so, the quality of the decoded image will suffer from the degradation. Considering this analogy, the use of decomposed loss might improve the decoded images for µ(x), encouraging the improvement of the orthonormality of the encoder/decoder in VAE.

heading: F.2 ADDITIONAL EXPERIMENTAL RESULT WITH OTHER CONDITION
text: In this Section, we provide the experimental results with other condition. We use essentially the same condition as described in Appendix D.2, except for the following conditions. The bottleneck size and λ are set to 256 and 10000, respectively. The encoder network is composed of CNN(9, 9, 2, 64, GDN) - CNN(5, 5, 2, 64, GDN) - CNN(5, 5, 2, 64, GDN) - CNN(5, 5, 2, 64, GDN) - FC(1024, 2048, softplus) - FC(2048, 256, None)×2 (for µ and σ) in encoder. The decoder network is composed of FC(256, 2048, softplus) - FC(2048, 1024, softplus) - CNN(5, 5, 2, 64, IGDN) - CNN(5, 5, 2, 64, IGDN) - CNN(5, 5, 2, 64, IGDN)-CNN(9, 9, 2, 3, IGDN).
Figures 24a and 24b show the averages of σj(x)−2 as well as the average and the standard deviation of 2βσj(x)
2D′j(z) in the conventional loss form and the decomposed loss form, respectively. When using the conventional loss form, the mean of 2βσj(x) 2D′j(z) is 1.25, which is closer to 1 than the
mean 1.83 in Section 4.2. This suggests that the implicit transform is closer to the orthonormal. The possible reason is that a bigger reconstruction error is likely to cause the interference to RD-trade off and a slight violation of the theory, and it might be compensated with a larger lambda. When using the decomposed loss form, the mean of 2βσj(x)
2D′j(z) is 0.95, meaning almost unit norm. These results also support that VAE provides the implicit orthonormal transform even if the lambda or bottleneck size is varied.

heading: G ADDITIONAL EXPERIMENTAL RESULT WITH MNIST DATASET
text: In this Appendix, we provide the experimental result of Section 4.2 with MNIST dataset3 consists of binary hand-written digits with a dimension of 768(=28× 28). We use standard training split which includes 50,000 data points. For the reconstruction loss, we use the binary cross entropy loss (BCE) for the Bernoulli distribution. We averaged BCE by the number of pixels.
The encoder network is composed of FC(768, 1024, relu) - FC(1024, 1024, relu) - FC(1024, bottleneck size) in encoder. The decoder network is composed of FC(bottleneck size, 1024, relu) - FC(1024, 1024, relu) - FC(1024, 768, sigmoid). The batch size is 256 and the training iteration number is 50,000. In this section, results with two parameters, (bottleneck size=32, λ=2000) and (bottleneck size=64, λ=10000) are provided. Note that since we averaged BCE loss by the number of pixels, β in the conventional β VAE is derived by 768/λ. Then, the model is optimized by Adam optimizer with the learning rate of 1e-3, using the conventional (not decomposed) loss form.
We use a PC with CPU Intel(R) Core(TM) i7-6850K CPU @ 3.60GHz, 12GB memory equipped with NVIDIA GeForce GTX 1080. The simulation time for each trial is about 10 minutes, including the statistics evaluation codes.
Figure 25 shows the averages of σj(x)−2 as well as the average and the standard deviation of 2 βσj(x) 2D′j(z). In both conditions, the means of 2 βσj(x)
2D′j(z) averages are also close to 1 except in the dimensions where σj(x)−2 is less than 10. These results suggest the theoretical property still holds when using the BCE loss. In the dimensions where σj(x)−2 is less than 10, the 2βσj(x) 2D′j(z) is somewhat lower than 1. The possible reason is that DKL(j)(·) in such dimension is 0 for some inputs and is larger than 0 in other inputs. The understanding of the transition region needs further study.
3http://yann.lecun.com/exdb/mnist/

heading: H DERIVATION/EXPLANATION IN RDO-RELATED EQUATION EXPANSIONS

heading: H.1 APPROXIMATION OF DISTORTION IN UNIFORM QUANTIZATION
text: Let T be a quantization step. Quantized values ẑj is derived as k T , where k = round(zj/T ). Then d, the distortion per channel, is approximated by
d = ∑ k ∫ (k+1/2)T (k−1/2)T p(zj)(zj − k T )2 dzj
' ∑ k T p(k T ) ∫ (k+1/2)T (k−1/2)T 1 T (zj − k T )2 dzj
= T 2
12 ∑ k T p(k T )
' T 2
12 . (67) Here, ∑ k T p(k T ) ' ∫∞ −∞ p(zj)dzj = 1 is used. The distortion for the given quantized value is
also estimated as T 2/12, because this value is approximated by ∫ (k+1/2)T (k−1/2)T 1 T (zj − k T ) 2 dzj .

heading: H.2 APPROXIMATION OF RECONSTRUCTION LOSS AS A QUADRATIC FORM.
text: In this appendix, the approximations of the reconstruction losses as a quadratic form tδxGxδx+Cx are explained for the sum of square error (SSE), binary cross entropy (BCE) and Structural Similarity (SSIM). Here, we have borrowed the derivation of BCE and SSIM from Kato et al. (2020), and add some explanation and clarification to them for convenience. We also describe the log-likelihood of the Gaussian distribution.
Let x̂ and x̂i be decoded sample Decθ(z) and its i-th dimensional component respectively. δx and δxi denote x− x̂ and xi− x̂i, respectively. It is also assumed that δx and δxi are infinitesimal. The details of the approximations are described as follows.

heading: Sum square error:
text: In the case of sum square error,Gx is equal to Im. This can be derived as:
m∑ i=1 (xi − x̂i)2 = m∑ i=1 δx2i = tδxImδx. (68)

heading: Binary cross entropy:
text: Binary cross entropy is a log likelihood of the Bernoulli distribution. The Bernoulli distribution is described as:
pθ(x|z) = m∏ i=1 x̂i xi (1− x̂i)(1−xi) (69)
Then, the binary cross-entropy (BCE) can be expanded as: − log pθ(x|z) = − log m∏ i=1 x̂i xi (1− x̂i)(1−xi)
= m∑ i=1 (−xi log x̂i − (1− xi) log (1− x̂i))
= ∑ i (−xi log(xi + δxi)− (1− xi) log(1− xi − δxi))
= ∑ i ( −xi log ( 1 + δxi xi ) − (1− xi) log ( 1− δxi 1− xi )) + ∑ i (−xi log(xi)− (1− xi) log(1− xi)). (70)
Figure 26: Graph of 12
( 1 x + 1 1−x ) in the BCE approximation.
Here, the second term of the last equation is a constant Cx depending on x. Using log(1 + x) = x− x2/2 +O(x3), the first term of the last equation is further expanded as follows:∑
i
( −xi ( δxi xi − δxi 2 2xi2 ) − (1− xi) ( − δxi 1− xi − δxi 2 2 (1− xi)2 ) +O ( δxi 3 ))
= ∑ i ( 1 2 ( 1 xi + 1 1− xi ) δxi 2 +O ( δxi 3 )) . (71)
As a result, a metric tensor Gx can be approximated as the following positive definite Hermitian matrix:
Gx =  1 2 ( 1 x1 + 11−x1 ) 0 . . . 0 12 ( 1 x2 + 11−x2 ) . . .
... ...
. . .  . (72)
Here, the loss function in each dimension 12 ( 1 x1 + 11−x1 ) is a downward-convex function as shown in Figure 26.

heading: Structural similarity (SSIM):
text: Structural similarity (SSIM) (Wang et al., 2001) is widely used for picture quality metric, which is close to human subjective quality. Let SSIM be a SSIM value between two pictures. The range of the SSIM value is between 0 and 1. The higher the value, the better the quality. In this appendix, we also show that (1− SSIM) can be approximated to a quadratic form such as tδx Gxδx. SSIMN×N(h,v)(x,y) denotes a SSIM value between N ×N windows in pictures X and Y , where x ∈ RN2 and y ∈ RN2 denote N × N pixels cropped from the top-left coordinate (h, v) in the images X and Y , respectively. Let µx, µy be the averages of all dimensional components in x, y, and σx , σy be the variances of all dimensional components in x, y in the N × N windows, respectively. Then, SSIMN×N(h,v)(x,y) is derived as
SSIMN×N(h,v)(x,y) = 2µxµy µx2 + µy2 · 2σxy σx2 + σy2 . (73)
In order to calculate a SSIM value for a picture, the window is shifted in a whole picture and all of SSIM values are averaged. Therefore, if ( 1− SSIMN×N(h,v)(x,y) ) is expressed as a quadratic form tδx G(h,v)x δx, (1− SSIM) can be also expressed in quadratic form tδx Gxδx.
Let δx be a minute displacement of x. µδx and σδx2 denote an average and variance of all dimensional components in δx, respectively. Then, SSIM between x and x + δx can be approximated as:
SSIMN×N(h,v)(x,x+ δx) ' 1− µδx
2
2µx2 − σδx
2
2σx2 +O
( (|δx|/|x|)3 ) . (74)
Then µδx2 and σδx2 can be expressed as
µδx 2 = tδxMδx, where M =
1
N2  1 1 . . . 1 1 1 . . . 1 ... ... . . .
... 1 1 . . . 1  , (75) and
σδx 2 = tδx V δx, where V =
1
N IN −M, (76) respectively. As a result, ( 1− SSIMN×N(h,v)(x,x+ δx) ) can be expressed in the following quadratic form as:
1− SSIMN×N(h,v)(x,x+ δx) ' tδx G(h,v)xδx, where G(h,v)x = ( 1
2µx2 M +
1
2σx2 V
) .
(77) It is noted thatM is a positive definite Hermitian matrix and V is a positive semidefinite Hermitian matrix. Therefore, G(h,v)x is a positive definite Hermitian matrix. As a result, (1 − SSIM) can be also expressed in quadratic form tδx Gxδx, whereGx is a positive definite Hermitian matrix.

heading: Log-likelihood of Gaussian distribution:
text: Gaussian distribution is described as:
pθ(x|z) = m∏ i=1 1√ 2πσ2 e−(xi−x̂i) 2/2σ2 = m∏ i=1 1√ 2πσ2 e−δxi 2/2σ2 , (78)
where σ2 is a variance as a hyper parameter. Then, the log-likelihood of the Gaussian distribution is denoted as:
− log pθ(x|z) = − log m∏ i=1 1√ 2πσ2 e−δx 2 i /2σ 2 = 1 2σ2 m∑ i=1 δx2i + m 2 log(2πσ2). (79)
The first term can be rewritten as (1/2σ2) tδxImδx. Thus,Gx = (1/2σ2) Im holds. Cx is derived as the second term of the last equation in Eq.78.

heading: H.3 DETAILED EXPLANATION OF KL DIVERGENCE AS A RATE OF ENTROPY CODING.
text: This appendix explains the detail how KL divergence can be interpreted as a rate in the transform coding. In the transform coding, input data is transformed by an orthonormal transform. Then, the transformed data is quantized, and an entropy code is assigned to the quantized symbol, such that the length of the entropy code is equivalent to the logarithm of the estimated symbol probability.
It is generally intractable to derive the rate and distortion of individual symbols in the ideal information coding. Thus, we first discuss the case of uniform quantization. Let Pzj and Rzj be the probability and rate in the uniform quantization coding of zj ∼ N (zj ; 0, 1). Here, µj(x) and σj(x)2 are regarded as a quantized value and a coding noise after the uniform quantization, respectively. Let T be a quantization step size. The coding noise after quantization is T 2/12 for the quantization step size T , as explained in Appendix H.1. Thus, T is derived as T = 2 √ 3σj(x) from σj(x)2 = T 2/12.
We also assume σj(x)2 1. As shown in Fig.27a, Pzj is denoted by ∫ µj(x)+T/2 µj(x)−T/2 p(zj)dzj where p(zj) is N (zj ; 0, 1). Using Simpson’s numerical integration method and ex = 1 + x + O(x2) expansion, Pzj is approximated as:
Pzj ' T
6
( p(µj(x) − T2 ) + 4p(µj(x)) + p(µj(x) + T 2 ) )
= Tp(µj(x))
6
( 4 + e 4µj(x)T−T 2 8 + e −4µj(x)T−T 2 8 ) ' Tp ( µj(x) ) ( 1− T 2/24
) = √ 6
π σj(x) e
−(µj(x)2)/2 (
1− σj(x)
2
2
) . (80)
Using log(1 + x) = x+O(x2) expansion, Rµσ is derived as:
Rzj = − logPzj ' 1
2
( µj(x) 2 + σj(x) 2 − log σj(x)2 − log 6
π
) = DKLj(x)(·) + 1
2 log
πe
6 . (81)
When Rzj and DKLj(x)(·) in Eq. 2 are compared, both equations are equivalent except a small constant difference 12 log(πe/6) ' 0.176 for each dimension. As a result, KL divergence for j-th dimension is equivalent to the rate for the uniform quantization coding, allowing a small constant difference.
To make theoretical analysis easier, we use the simpler approximation as Pzj = T p(µj(x)) = 2 √
3σj(x) p(µj(x)) instead of Eq.80, as shown in Fig.27b. Then, Rzj is derived as:
Rzj = − log(2 √ 3 σj(x) p(µj(x))) = Eq. 6 + 1
2 log
πe
6 . (82)
This equation also means that the approximation of KL divergence in Eq. 6 is equivalent to the rate in the uniform quantization coding with Pzj = 2 √ 3σj(x) p(µj(x)) approximation, allowing the same small constant difference as in Eq. 81. It is noted that the approximation Pzj = 2 √
3σj(x) p(µj(x)) in Figure 27b can be applied to any kinds of prior PDFs because there is no explicit assumption for the prior PDF. This implies that the theoretical discussion after Eq. 6 in the main text will hold in arbitrary prior PDFs.
Finally, the meaning of the small constant difference 12 log πe 6 in Eqs. 81 and 82 is shown. Pearlman & Said (2011) explains that the difference of the rate between the ideal information coding and uniform quantization is 12 log πe 6 . This is caused by the entropy difference of the noise distributions. In the ideal case, the noise distribution is known as a Gaussian. In the case the noise variance is σ2, the entropy of the Gaussian noise is 12 log(σ
22πe). For the uniform quantization with a uniform noise distribution, the entropy is 12 log(σ 212). As a result, the difference is just 12 log πe 6 . Because the rate estimation in this appendix uses a uniform quantization, the small offset 12 log πe 6 can be regarded as a difference between the ideal information coding and the uniform quantization. As a result, KL divergence in Eq. 2 and Eq. 6 can be regarded as a rate in the ideal informaton coding for the symbol with the mean µj(x) and variance σj(x)2.

Review summaries:
The paper builds on a branch of recent works that consider and analyse Variational autoencoders (VAE) from the view point of data compression. This started from Alemi et al. 2018, where the authors consider and analyse the mutual information between data and latent codes and culminates in Kato et al. 2020, where the authors consider models in which both the encoder and the decoder are assumed as deterministic (isometric) mappings. The submitted paper aims at reconciling this type of models with standard VAEs by claiming to prove that VAEs can be obtained from the former by an non-linear component-wise scaling of the latent space. The authors claim among others that this approach allows to estimate the data probability density from the learned VAE model.