 # Maximum likelihood estimation of the Fisher-Bingham distribution via efficient calculation of its normalizing constant

This paper proposes an efficient numerical integration formula to compute the normalizing constant of Fisher–Bingham distributions. This formula uses a numerical integration formula with the continuous Euler transform to a Fourier-type integral representation of the normalizing constant. As this method is fast and accurate, it can be applied to the calculation of the normalizing constant of high-dimensional Fisher–Bingham distributions. More precisely, the error decays exponentially with an increase in the integration points, and the computation cost increases linearly with the dimensions. In addition, this formula is useful for calculating the gradient and Hessian matrix of the normalizing constant. Therefore, we apply this formula to efficiently calculate the maximum likelihood estimation (MLE) of high-dimensional data. Finally, we apply the MLE to the hyperspherical variational auto-encoder (S-VAE), a deep-learning-based generative model that restricts the latent space to a unit hypersphere. We use the S-VAE trained with images of handwritten numbers to estimate the distributions of each label. This application is useful for adding new labels to the models.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

### 1.1 Fisher–Bingham distribution

The Fisher–Bingham distribution is defined as a multivariate normal distribution restricted on a unit sphere.

###### Definition 1

For a -dimensional multivariate normal distribution with a mean

and a variance-covariance matrix

, the Fisher–Bingham distribution is given by the density function

 f(x;μ,Σ):=1Cexp(−xTΣ−1x2+xTΣ−1μ)dSp−1(x),

where and

 C=C(Σ−12,Σ−1μ) := ∫Sp−1exp(−xTΣ−1x2+xTΣ−1μ)dSp−1(x)

is the normalizing constant and is the uniform measure in the -dimensional sphere .

The Fisher–Bingham distribution plays an essential role in directional statistics, which is concerned with data on various manifolds, especially data represented in a high-dimensional sphere. For example, wind direction and the geomagnetic field are common types of data that can be represented on a sphere . In addition, data on a hypersphere are used in link prediction of networks and image generation. Therefore, the Fisher–Bingham distribution, a normal distribution restricted on a unit sphere, is commonly used in this field.

However, the spherical domain causes some problems when using Fisher–Bingham distributions. One such problems is calculating the normalizing constant. As it is difficult to calculate it analytically, a numerical method is necessary. The saddlepoint approximation method is a numerical method for computing the normalizing constant developed by kume2005saddlepoint. Another approach, the holonomic gradient method considered by kume2018exact, computes the normalizing constant as well. However, these methods have some limitations. The saddlepoint approximation method is not as accurate as the holonomic gradient method, which is theoretically exact because the problem of calculating is mathematically characterized by solving an ODE. However, the holonomic gradient method is computationally expensive and cannot be applied to calculate the normalizing constant of high-dimensional distributions. Hence, it is necessary to create a numerical method that is efficient, numerically stable, and accurate.

To construct such a numerical method, the following details about Fisher–Bingham distributions are required (kume2018exact).

Since any orthogonal transformation in is isometric, the parameter dimensions are reduced from to

by singular value decomposition. Therefore, we have

 C(Σ−12,Σ−1μ)=C(Δ−12,Δ−1Oμ),

where and

is the orthogonal matrix obtained from

. Thus, without loss of generality, we can assume that the variance-covariance matrix is diagonal. After reducing the parameter dimensions to , the normalizing constant becomes

 C(Δ−12,Δ−1Oμ)=C(θ,γ) := ∫Sp−1exp(p∑i=1(−θix2i+γixi))dSp−1(x),

where

 θ=(θ1,⋯,θp)=(12δ21,⋯,12δ2p)=diag(Δ−12)

and

 γ=(γ1,⋯,γp)=Δ−1Oμ.

Since is restricted on a unit sphere, we have

 C(θ+cI,γ) = ∫Sp−1exp(p∑i=1(−(θi+c)x2i+γixi))dSp−1(x) = ∫Sp−1exp(−c+(p∑i=1(−θix2i+γixi)))dSp−1(x) = e−c∫Sp−1exp(p∑i=1(−θix2i+γixi))dSp−1(x) = e−cC(θ,γ),

where is a real number and . If we put

 f(x;θ,γ):=1C(θ,γ)exp(p∑i=1(−θix2i+γixi))dSp−1(x).

then we have

 f(x;θ+cI,γ) = 1C(θ+cI,γ)exp(p∑i=1(−(θi+c)x2i+γixi))dSp−1(x) = ecC(θ,γ)exp(−c+(p∑i=1(−θix2i+γixi)))dSp−1(x) = 1C(θ,γ)exp(p∑i=1(−θix2i+γixi))dSp−1(x) = f(x;θ,γ).

As a result, if the normalizing constant is obtained, can also be obtained. Moreover, for the maximum likelihood estimation (MLE), as , can be shifted to for all . Additionally, because the unit sphere is symmetrical,

 C(θ,|γ|) =∫Sp−1exp(p∑i=1(−θix2i+|γi|xi))dSp−1(x) =∫Sp−1exp(p∑i=1(−θix2i+γixi))dSp−1(x) =C(θ,γ).

As a result, it can be assumed that has non-negative entries when calculating the normalizing constant.

### 1.2 Aim of this paper

In this paper,

1. we propose an efficient numerical integration formula to compute the normalizing constant.

2. we apply this formula to perform MLE.

3. we apply MLE to the latent variables of the hyperspherical variational auto-encoder (S-VAE) (davidson2018hyperspherical).

The normalizing constant of Fisher–Bingham distributions can be represented in a Fourier integration form. Therefore, we can use the numerical integration formula with the continuous Euler transform introduced by ooura2001continuous. Note that the continuous Euler transform is useful for calculating the normalizing constant and MLE.

This method can be applied to the MLE of high-dimensional data, such as the latent variables of S-VAE (davidson2018hyperspherical

), a generating model used in machine learning. The dimensions of the hyperspherical variational auto-encoder rely on the complexity of the data. For example, for human face data, there may be 100 dimensions of the latent variables.

### 1.3 Organization of this paper

This paper is organized as follows. In Section 2, we make some general remarks about the Fisher–Bingham distribution and the Fourier transform representation of the normalizing constant. In Section 3, we explain the continuous Euler transform and its use for numerical computation of the normalizing constant. In Section 4, we discuss the calculation of the gradient of the normalizing constant, which is necessary for MLE. Subsequently, the MLE algorithm is provided. In Section 5, we demonstrate some MLE numerical experiment to show the effectiveness of this method. In Section 6, we show the application of MLE in the S-VAE whose latent space includes high-dimensional data on a hypersphere.

## 2 Fourier transform representation of the normalizing constant

### 2.1 Laplace inversion representation

This section explains how the normalizing constant can be represented in a simpler form, as derived by kume2005saddlepoint and kume2018exact. The first step is to change the range of integration of the normalizing constant from a -dimensional hypersphere to a one-dimensional line . The integration range is then shifted to . Note that the derivation can be somewhat technical, but the calculation itself is not difficult.

First, the distribution of

independent normal random variables

is

 f(x1,⋯,xp)=∏pi=1θ12iπp2exp(−p∑i=1θi(xi−μi)2).

We then apply the variable transform

 ⎧⎪ ⎪⎨⎪ ⎪⎩r=∑pi=1x2i=xTxϕ=(ϕ1,⋯,ϕp)=(x1r12,⋯,xpr12)=xr12

to and integrate it with respect to . Then, the marginalized distribution becomes

 fmrg(r)= 12π−p2(p∏i=1θ12i)^C(rθ,r12γ) ×exp(−14p∑i=1γ2iθi)rp2−1, (1)

where

 γ=(2θ1μ1,⋯,2θpμp)

and

 ^C(rθ,r12γ) = ∫Sp−1exp(−p∑i=1(rθiϕ2i−r12γiϕi))dSp−1(ϕ). (2)

When , Equation (2) matches the definition of the normalizing constant . As a result, based on Equation (1), we obtain

 C(θ,γ)=2πp2(p∏i=1θ−12i)fmrg(1)exp(14p∑i=1γ2iθi). (3)

Therefore, if the distribution can be represented in a one-dimensional integration form, the goal will be achieved.

is

 L(t)=exp(∑pi=1(γ2i4(θi+t)−γ2i4θi))∏pi=1√1+tθi.

Since , the moment generating function is the same as the Laplace transform of . Thus, with the inverse Laplace transform, we obtain

 fmrg(r)=12πi∫iR+t0L(t)ertdt, (4)

where . Substituting Equation (4) into Equation (3), we get

 C(θ,γ)=−iπp2−1∫iR+t0p∏i=1exp(γ2i4(θi+t))√θi+tetdt. (5)

This is the Laplace inversion representation of the normalizing constant, which contains a one-parameter integration of . The Fourier transform representation can be easily deduced by applying the variable transform to Equation (5).

### 2.2 Fourier-type integral representation

###### Theorem 1
 C(θ,γ)=πp2−1e−t0∫RA(t;θ,γ)eitdt, (6)

where

 A(t;θ,γ)=p∏i=1exp(γ2i4(θi−it−t0))√θi−it−t0.

Now, we derive the representation of the normalizing constant in the form of a Fourier integral. Therefore, the next step is to apply a numerical integration formula to calculate this Fourier integral.

## 3 Continuous Euler transform for numerical calculation

### 3.1 Continuous Euler transform

Since the normalizing constant is represented as a Fourier-type integral in Equation (6), it is necessary to find an accurate numerical integration formula to calculate this integration. Here, note that the function decays slowly. With a usual trapezoidal formula, the slower the function decays, the slower the convergence of the numerical integration becomes. It is necessary to reinforce the decay of the function to improve the accuracy. Here, we can adopt the continuous Euler transform by ooura2001continuous, which can accelerate the convergence of the Fourier transform. Moreover, tanaka2014error proved that adding the continuous Euler transformation to the trapezoidal formula made the integration converge rapidly. The details of the continuous Euler transform have been provided by ooura2001continuous and tanaka2014error.

We apply the continuous Euler transform to Equation (6). Choose and to satisfy and . Choose . Let be an integer with

 N≥2d(ωd+ωu)ω2uπω2d.

Let , and be defined by

 h= ⎷2πd(ωd+ωu)ω2dN,p=√Nhωd,q=√ωdNh4.

Then, we get

 C(θ,γ) = πp2−1e−t0∫RA(t;θ,γ)eitdt ≈ πp2−1e−t0∫Rw(|t|,p,q)A(t;θ,γ)eitdt (7) ≈ πp2−1e−t0hN∑n=−N−1w(|nh|;p,q)A(nh,θ,γ)einh (8) =: C(N,h)w(θ,γ), (9)

where

 w(x;p,q)=12erfc(xp−q).

Equation (8) uses the trapezoidal formula.

The accuracy of the numerical formula is prooved using Theorem 2.

###### Theorem 2
 |C(θ,γ)−C(N,h)w(θ,γ)|≤Poly(N)exp⎛⎜⎝− ⎷πdω2dN2(ωd+ωu)⎞⎟⎠,

where is a polynomial of .

With this theorem, if the normalizing constant is approximated by , the error converges to as . In other words, any accuracy can be achieved if a practically large enough is taken. For instance, if the normalizing constant with an error less than is necessary, is sufficient. Subsequently, if the parameter dimensions are about 100, it only takes 20 ms, whereas the holonomic gradient method takes about 15 s with 10 dimensions (sei2015calculating). The details of the numerical experiments are provided in Section 5. An efficient numerical integration formula for the normalizing constant can be obtained using this method.

## 4 MLE optimization using the continuous Euler transform

The likelihood function of the Fisher–Bingham distribution has been provided by kume2018exact. With a observed data matrix , put and , the likelihood function is

 logL(Σ−12,Σ−1μ,X) = logn∏i=1⎛⎜ ⎜ ⎜ ⎜⎝exp(xTiΣ−1μ−xTiΣ−1xi2)C(Σ−12,Σ−1μ)⎞⎟ ⎟ ⎟ ⎟⎠ = −nlogC(Δ−12,γ)−n∑i=1(xTiΣ−12xi−xTiΣ−1μ) = −nlogC(Δ−12,γ)−n%tr(AOTΔ−12O−OBγT) = −n(logC(θ,γ)+tr(AOTdiag% (θ)O+OBγT))

with

 Σ−1=OTΔ−1O, C(Σ−12,Σ−1μ)=C(Δ−12,γ)=C(θ,γ), γ=Δ−1Oμ

and

 Δ−12=diag(θ).

Therefore, maximizing the likelihood function

 logL(Σ−12,Σ−1μ,X)

is equivalent to minimizing

 logL(θ,γ,O) := logC(θ,γ)+tr(AOTdiag(θ)O+OBγT). (10)

In this section, is also called the likelihood function, although

 logL(θ,γ,O)=−1nlogL(Σ−12,Σ−1μ,X).

Thus, it is possible to optimize Equation (10) by iteratively updating the parameters that decrease the likelihood value. First, we consider the optimization problem in and for a fixed . Then, the optimal can be obtained by minimizing mantaining the value of and fixed,. The formula for updating has been published by kume2018exact. Although the outline for updating is mentioned in this section,

is fixed as an identity matrix, that is,

is assumed to be diagonal in numerical experiments.

It is necessary to calculate the partial derivatives of when optimizing and .

 ∂logL(θ,γ,O)∂θ=∂C(θ,γ)∂θ1C(θ,γ)+diag(OAOT) (11)

and

 ∂logL(θ,γ,O)∂γ=∂C(θ,γ)∂γ1C(θ,γ)+BTOT. (12)

The partial derivatives of are needed to calculate Equation (11) and Equation (12).

In Theorem 1, the Fourier transform representation of the normalizing constant becomes

 C(θ,γ)=πp2−1et0∫RA(t;θ,γ)eitdt,

where

 A(t;θ,γ)=p∏i=1exp(γ2i4(θi−it−t0))√θi−it−t0.

By changing the order of integration and differentiation, the partial derivatives of become

 Cθi(θ,γ) :=∂C(θ,γ)∂θi =πp2−1et0∫R∂A(t;θ,γ)∂θieitdt, Cγi(θ,γ) :=∂C(θ,γ)∂γi =πp2−1et0∫R∂A(t;θ,γ)∂γieitdt.

The continuous Euler transform can also be applied to these calculations to numerically calculate the derivatives efficiently.

If we put

 = {exp(γ2i4(θi−it−t0))12(θi−it−t0) ×(−γ2i2(θi−it−t0)−1)}∏j≠iexp(γ2i4(θi−it−t0))√θi−it−t0, = exp(γ2i4(θi−it−t0))γi2(θi−it−t0)32∏j≠iexp(γ2i4(θi−it−t0))√θi−it−t0.

Then, we get

 Cθi(θ,γ) = πp2−1et0∫RAθi(t;θ,γ)eitdt ≈ (13)
 Cγi(θ,γ) = πp2−1et0∫RAγi(t;θ,γ)eitdt ≈ (14)

Therefore, the derivatives of the normalizing constant, as well as the derivatives of the likelihood function, can be calculated using Equation (11) and (12). As a result, it is possible to optimize and with a fixed .

Given that the algorithm to optimize with fixed and values was developed by kume2018exact, this will not be demonstrated in this paper.

 A=diag(θ)OAOT−OAOTdiag(θ)+γBTOT (15)

and

 ^v=A−AT. (16)

must be symmetrical for to be an optional orthogonal matrix (kume2018exact). Moreover, if is symmetrical, a curve with reduces to a single point; this can be used as a stopping criterion. The proof of this stopping criterion was provided by kume2018exact.

Now, we have obtained all parts necessary to optimize the likelihood function. Consequently, we can discuss the approach for obtaining the MLE of Fisher–Bingham distributions. This algorithm performs the same steps as the algorithm given by kume2018exact, but the method used to calculate the derivatives of the likelihood function and the normalizing constant are different. Although Kume and Sei used the holonomic gradient method, the continuous Euler transform is adopted in this paper.

###### Algorithm

Update the given , , and as follows until the differentiation of , , and becomes small enough:

1.  ^θ=θ+∂logL(θ,γ,O)∂θδθ.

The partial derivatives of the likelihood function are obtained by substituting Equation (13) into Equation (11), and is a real number such that

 logL(^θ,γ,O)
2.  ^γ=γ+∂logL(θ,γ,O)∂γδγ.

The partial derivatives of the likelihood function are obtained by substituting Equation (14) into Equation (12), and is a real number such that

 logL(θ,^γ,O)
3.  ^O=e^vt0O.

is obtained by substituting Equation (15) into Equation (16), and is a real number such that

 logL(θ,γ,^O)

This algorithm is based on the gradient descent method that converges linearly. If faster convergence is required, it is preferable to use the quasi-Newton method.

## 5 Numerical experiments

### 5.1 Normalizing constant

We compare the calculation method for the normalizing constant of the Bingham distribution, the holonomic gradient method, and the saddlepoint approximation method with our method. The Bingham distribution is a special case of the Fisher–Bingham distribution, which fixes .

###### Definition 2

For a -dimensional multivariate distribution, the Bingham distribution is given by the density function

 f(x;θ):=1C(θ)exp(−p∑i=1θix2i)dSp−1(x),

where and

 C(θ):=∫Sp−1exp(−p∑i=1θix2i)dSp−1(x)

is the normalizing constant and is the uniform measure in the -dimensional sphere .

Theorem 2 is used to compute the normalizing constant of the Bingham distribution. More precisely, the numerical integration formula is

 C(N,h)w(θ)=πp2−1et0hN∑n=−N−1w(|nh|;p,q)A(nh,θ)einh,

where

 A(t;θ)=p∏i=11√θi−it−t0.

The same parameters were used by sei2015calculating. The results of the holonomic gradient method and the saddlepoint approximation method were obtained from sei2015calculating. In our method, the parameter is fixed to 200. The holonomic gradient method must be theoretically exact because the problem is mathematically characterized by solving an ODE with numerical methods.

The accuracy of the holonomic gradient method and our method is very high. As shown in Table 1 to Table 4, the normalizing constants with different parameters calculated by these methods (the hg and ce columns) coincide until they reach 6 decimal places. However, the saddlepoint approximation method is not as good because the error is larger than . Besides, the calculation time for our method with 100 dimensions is about 20 ms, whereas the holonomic gradient method takes about 15 s with 10 dimensions (sei2015calculating).

In Table 1, columns 2 and 3 compare the saddlepoint approximation (spa) and the holonomic gradient method (hg) of with our method (ce); columns 4 and 5 compare the same quantities for .

The complex Bingham distribution, which is defined on a unit complex sphere, can be calculated analytically (kume2005saddlepoint).

###### Definition 3

For a -dimensional multivariate distribution, the complex Bingham distribution is given by the density function

 f(x;θ):=1C(θ)exp(−p∑i=1θix2i)dSp−1(x),

where and

 C(θ):=∫Sp−1exp(−p∑i=1θix2i)dSp−1(x)

is the normalizing constant and is the uniform measure in the -dimensional complex sphere .

In Table 2, we compare the saddlepoint, exact, hg, and our method for the complex Bingham distribution with parameters , that is, .

In Table 3, columns 2 and 3 compare the saddlepoint approximation (spa) and the holonomic gradient method (hg) of with our method (ce); columns 4 and 5 compare the same quantities for .

In Table 4, we compare the saddlepoint, exact, hg, and our method for the complex Bingham distribution with parameters , that is, .

The advantage of our method is its efficiency in high-dimensional cases. The calculation time of the normalizing constant with multiple dimensions is demonstrated in Figure 1. The parameters and are determined by the random number generator in C++ library and .

The calculation is very rapid and stable. The calculation time increases linearly. with the number of dimensions. Since the derivatives of the normalizing constant are calculated with the same numerical integration formula, we prove that the derivatives can be calculated efficiently. The MLE algorithm, introduced in Section 4, is based on the gradient descent method. The bottleneck depends on the optimizing method because the calculation of the gradients is not expensive. Since the gradient descent method converges linearly, it is not efficient enough for high-dimensional data. This problem can be solved using the quasi-Newton method.

### 5.2 Mle

In this section, some numerical experiments on the MLE using data with multiple dimensions are shown. The data are obtained by rejection sampling, where is assumed to be diagonal, that is, is not estimated.

#### 5.2.1 2-dimensional data

For data matrix (Figure 2), where

 xi iid∼f(x;θ∗,γ∗)=1C(θ∗,γ∗)exp(2∑i=1(−θ∗ix2i+γ∗ixi)), θ∗ =(1.02.0)andγ∗=(1.02.0),

the following MLE values

 ^θ =(1.045191.95468)and^γ=(1.046432.01098)

are obtained. Figure 2: Histogram of the sampling points on a unit circle

#### 5.2.2 3-dimensional data

For a data matrix (Figure 3), where

 xi iid∼f(x;θ∗,γ∗)=1C(θ∗,γ∗)exp(3∑i=1(−θ∗ix2i+γ∗ixi)), θ∗ =⎛⎜⎝1.02.03.0⎞⎟⎠andγ∗=⎛⎜⎝1.02.03.0⎞⎟⎠,

the following MLE values

 ^θ =⎛⎜⎝1.022151.996362.98587⎞⎟⎠and^γ=⎛⎜⎝0.904622.017932.99908⎞⎟⎠

are obtained.

#### 5.2.3 10-dimensional data

For a data matrix , where

 xiiid∼ f(x;θ∗,γ∗) = 1C(θ∗,γ∗)exp(10∑i=1(−θ∗ix2i+γ∗ixi)), θ∗ =(1,2,3,4,5,6,7,8,9,10)T, γ∗ =(1,2,3,4,5,6,7,8,9,10)T

the following MLE values

 ^θ=( 1.2849,2.6223,3.0729,4.3598,5.0906, 5.6851,6.2268,7.2384,7.9931,9.0433)T, ^γ=( 0.9962,2,1081,2,9225,4.1555,4.8213, 5.6873,6.3327,7.4991,8.1197,9.2781)T

are obtained.

## 6 Application of Fisher–Bingham distribution to S-VAE

### 6.1 Introduction to S-VAE

As mentioned in Section 5, our method enables to perform the MLE with high-dimensional data which cannot be achieved using other methods such as the holonomic gradient method. Because of this advantage, we can apply our method to the MLE of the latent variables of variational auto-encoders (VAE).

It is essential to comprehend the auto-encoders (AE) mechanism before introducing the VAE. AE is a generating model with two networks. One is called the encoder, while the other is called the decoder. The encoder transforms images into vectors with smaller dimensions than those of the images. For example, MNIST is a dataset with handwritten numbers from 0 to 9, and each image has 28

28 dimensions. With the encoder of the AE, an image with 28 28 dimensions is transformed into a vector with 2-20 dimensions. If an image becomes complicated, the dimensions of the latent space increase. For instance, the latent space of human face data has about 100 dimensions.

VAE is a generating model derived from AE. The difference between the VAE and AE is that the VAE assumes that the vectors obtained from images are generated from some distributions (doersch2016tutorial; kingma2013auto

). The encoder estimates the parameters of distributions, while the decoder generates an image from the vector, which is a sample of the estimated distributions. For example, in most cases, the vectors are assumed to be generated from some normal distributions in Euclidean space. In this section, the VAE implies that the normal distribution in Euclidean space is assumed. Considering MNIST, 10 normal distributions match with different numbers from 0 to 9. If we decode a vector obtained from the distribution of zero, something like zero will probably be generated. Therefore, the encoder of the VAE transforms images into parameters, such as the means and variations of the distributions, instead of vectors themselves. The advantage of the VAE compared to the AE the VAE that it gives the structure of the latent variables.

S-VAE is a VAE assuming that the vectors are generated from distributions restricted to a high-dimensional unit hypersphere. Moreover, the vectors are assumed to be generated from the von Mises–Fisher distribution.

###### Definition 4

For a -dimensional multivariate distribution, the von Mises–Fisher distribution is given by the density function

 f(x;μ,κ):=1C(κ)exp(κμTx)dSp−1(x),

where and

 C(κ):=∫Sp−1exp(κμTx)dSp−1(x)

is the normalizing constant and is the uniform measure in the -dimensional sphere .

The von Mises–Fisher distribution is a special case of the Fisher–Bingham distribution when restricting for all . When

, the von Mises–Fisher distribution becomes a uniform distribution on a unit sphere.

davidson2018hyperspherical

proposed two reasons to assume that the latent space is a unit sphere rather than Euclidean space. First, S-VAE allows for a uniform distribution on the hypersphere to be a prior, which is a truly uninformative prior. Since uniform distribution in Euclidean space does not exist, the VAE must make some informative assumptions about the prior. Second, with high-dimensional data, ”the soap-bubble-effect”, where the latent variables converge on a hyperspherical shell, is observed. Therefore, the assumption of the latent space as a Euclidean hyperplane would be improper because

is not homeomorphic to .

The AE, VAE and S-VAE are unsupervised machine learning, which means that the labels of images are not used for the training of networks. However, images with high similarities, such as several handwritten images of the number one, become close to each other in latent space. Therefore, some clusters are generated in latent space as shown in Figures 6 and 7 (davidson2018hyperspherical). Each color matches with each label from 0 to 9. Figure 7 is the Hammer projection, a projection that shows a sphere on a 2-dimensional plain. Therefore, the latent variable representation of each image can be obtained with a pre-trained model. The parameters of these clusters can then be estimated using MLE. The algorithm mentioned in Section 4 can be applied specifically to S-VAE.

In addition, the conditional VAE (CVAE) is a supervised model that uses the label information to train the networks.

It is common to consider adding more details to the labels of CVAE. For example, if a pre-trained CVAE with human faces is labeled with continents, such as Asia, Europe and Africa, we want to add the details of continents to the labels. Additionally, we must consider the situation in which a pre-trained CVAE with images of fish labeled with different species. When a new species is discovered, we may want to add this new species label to the network. In most cases, it is common to train the networks with the dataset to which new labels are appended. However, the computing cost is high. For example, the training of the VAE of MNIST with the latent space of a using a CPU (Dell latitude 5480) takes several hours in our experiment.

If the latent variables are assumed to be generated with Fisher–Bingham distributions, the distribution of the images with new labels can be obtained by estimating its parameters via MLE. In this way, since the computation cost of MLE is much lower than the training process of CVAE, new labels can be appended much more easily. In the VAE of MNIST, the estimation of the parameters takes only seconds. Therefore, the MLE algorithm mentioned in Section 4 can be applied for S-VAE. Attention should be paid to other methods for the MLE of the Fisher–Bingham distribution, such as the holonomic gradient method and the saddlepoint approximation method, are not appropriate for high-dimensional data.

### 6.2 Application results

In this section, the handwritten images dataset MNIST is used to train S-VAE. The experiments with the latent space and are shown. The models are trained according to davidson2018hyperspherical

. It is necessary to add a sigmoid layer to the decoder. The training process is increased from 1000 to 2000 epochs with early-stopping at 50. Besides these changes, the training procedures are mostly the same as those mentioned by

davidson2018hyperspherical. Note that the experiment results are not accurate enough. In this section, the accuracy of the model means roughly the degree of the similarity of generated images to the handwritten numbers. The similarity is judged by the recognition of the authors. However, some concerns remain, such as the accuracy of the models and that of the Fisher–Bingham distribution assumption.

#### 6.2.1 Latent space S2

The images generated by S-VAE are shown in Figure 9.