Nonparametric Bayesian Deep Visualization
Haruya Ishizuka
12

1

and Daichi Mochihashi2

Bridgeston Corporation haruya.ishizuka@bridgestone.com
The Institute of Statistical Mathematics daichi@ism.ac.jp

Abstract. Visualization methods such as t-SNE [1] have helped in knowledge discovery from high-dimensional data; however, their performance
may degrade when the intrinsic structure of observations is in low-dimensional
space, and they cannot estimate clusters that are often useful to understand the internal structure of a dataset. A solution is to visualize thelatent coordinates and clusters estimated using a neural clustering model.
However, they require a long computational time since they have numerous weights to train and must tune the layer width, the number of latentdimensions and clusters to appropriately model the latent space. Additionally, the estimated coordinates may not be suitable for visualization
since such a model and visualization method are applied independently.
We utilize neural network Gaussian processes (NNGP) [2] equivalent toa neural network whose weights are marginalized to eliminate the necessity to optimize weights and layer widths. Additionally, to determine
latent dimensions and the number of clusters without tuning, we propose a latent variable model that combines NNGP with automatic relevance determination [3] to extract necessary dimensions of latent space
and infinite Gaussian mixture model [4] to infer the number of clusters.
We integrate this model and visualization method into nonparametric
Bayesian deep visualization (NPDV) that learns latent and visual coordinates jointly to render latent coordinates optimal for visualization.
Experimental results on images and document datasets show that NPDV
shows superior accuracy to existing methods, and it requires less training
time than the neural clustering model because of its lower tuning cost.
Furthermore, NPDV can reveal plausible latent clusters without labels.
Keywords: Data visualization · Gaussian processes · Nonparametric
Bayesian models · Neural network

1

Introduction

Visualization methods such as t-SNE [1], which compress input to two- or threedimensional visual coordinates to be mapped on a scatter plot, provide a useful
overview of high-dimensional data, and a number of methods have been proposed. These methods estimate visual coordinates based on the similarity ofdata points and fall into two categories. The first category is local methods topreserve neighbor structures in the original space [5–9]. These methods utilize
nearest neighbor graph that represents pairwise similarity to retain distances between neighbors. Among them, t-SNE [1] and UMAP [10] are the most popular2

H. Ishizuka and D. Mochihashi

Fig. 1: 3D Visualization of 100-dimensional data generated by transforming the
mammoth data in (a). Since NPDV(MF) estimates not only visual coordinates
but also latent clusters, the associated plot in (b) can be colored by the cluster
assignments, unlike existing methods in (c)–(e). See section 5 for details.
algorithms. The second is global methods to exploit relationships among threepoints, and preserve distances between data points distant from one another [11–13]. Particularly, Trimap [14] arguably shows comparable accuracy to t-SNE and
UMAP. There is a trade-off between preserving the local or global structures.
PaCMAP [15] exceptionally preserves both by combining the loss functions oflocal and global methods. They have improved knowledge discovery in various
domains, such as bioinformatics [16] and audio processing [17].
However, they have several drawbacks to reveal hidden structures behinddatasets. First, their performance may degrade when the observations are distributed on a low-dimensional manifold embedded in the observation space. In
this case, the similarity between observations may differ from that in the manifold, resulting in an inaccurate visualization. This problem worsens when the
manifold is embedded by a highly nonlinear function. Second, they cannot estimate clusters. Visualization together with clusters provides an intuitive understanding of the internal structure of a dataset; however, most of these methods
only estimate visual coordinates. Supervised dimensionality-reduction [18–20]
utilize labels to address this drawback; however, these are not always available.
We present these problems through a simulation example. In this example, the3D mammoth data shown in Fig.1 (a) is embedded into a 100-dimensional space
nonlinearly by a neural network; then, this data is visualized in 3D space using
several methods3 . Exiting methods in Fig.1 (c)–(e) evidently fail to recover theoriginal mammoth shape. In addition, the lack of estimating clusters make it
difficult to interpret the internal structure of resulting plot.
A solution here is to visualize latent coordinates and clusters estimated using a neural clustering model [21–23]; however, this approach presents otherissues. As these models can accurately model a low-dimensional manifold byleveraging a neural autoencoder and perform clustering simultaneously, this approach addresses the aforementioned issues. However, they often require a longcomputational time to optimize model performance since they have numerous
neural weights to train and need to search their appropriate hyperparametersettings, which is not suitable for scientific visualization. Particularly, the widthof hidden layers, the number of latent dimensions, and clusters that are criti3

Rotatable plots are provided as an html file in the supplemental material.

Nonparametric Bayesian Deep Visualization

3

Fig. 2: Visualization flow of NPDV. NPDV integrates the estimation of latentcoordinates X from observations Y using NN-iWMM and that of visual coordinates V from X. Clustering is also performed by NN-iWMM.
cal hyperparameters to reveal the latent structure. Furthermore, the estimatedlatent coordinates may not be suitable for visualization as such a model andvisualization method are independently applied.
We utilize neural network Gaussian processes (NNGP) [2], which are equivalent to a neural network whose weights are marginalized, to address these issues.
The marginalization greatly reduces the computational time by eliminating thenecessity to optimize numerous weights and layer width while exploiting thepower of neural networks. Additionally, to determine the number of latent dimensions and clusters without tuning, we propose the neural network infinite
warped mixture model (NN-iWMM) by combining NNGP with automatic relevance determination [3] to extract the necessary dimensions of latent space andthe infinite Gaussian mixture model [4] to infer the number of clusters. Finally,
NN-iWMM and visualization methods are integrated as shown in Fig. 2, into
nonparametric Bayesian deep visualization (NPDV) that jointly infers latentand visual coordinates to render latent coordinates optimal for visualization.
Based on NPDV, we introduce NPDV(MF), which employs matrix factorization to linearly reduce the dimensionality, and NPDV(t-SNE) based on t-SNE;
both methods enables to visualize the internal structure of dataset by utilizing
the estimated clusters. As shown in Fig.1 (b), we can observe from the simulation study that NPDV(MF) could accurately recover the mammoth shape.
NPDV(t-SNE) achieves better accuracy than existing methods and NPDV(MF)for real-world data. Furthermore, NPDV(t-SNE) shows two preferable properties in unsupervised settings: (1) it takes considerably less training time than
the neural clustering model and (2) has the ability to reveal plausible clusters
without label information..
The remainder of this paper is organized as follows. We introduce the preliminaries in Section 2, and the proposed models and their training algorithm in
Sections 3 and 4, respectively. Subsequently, we demonstrate their advantages
through simulation and real data experiments in Sections 5 and 6.

2

Infinite Warped Mixture Model

The proposed models are based on the infinite warped mixture model (iWMM)[24], an extension of the Gaussian process latent variable model (GPLVM) [25]

4

H. Ishizuka and D. Mochihashi

that uses Gaussian processes [26] for dimensionality reduction. We introduce thenotation and iWMM as preliminary. Let D, Q, and S denote the dimensionalities
Q N
of observations Y = {yi ∈ RD }N
i=1 , latent coordinates X = {xi ∈ R }i=1 , andS Nvisual coordinates V = {vi ∈ R }i=1 , respectively. Q is smaller than D andlarger than S, and S is typically set to two or three. N represents the numberof observations. N (m, C) is a multivariate Gaussian distribution with a meanm and covariance C. IN represents an N -dimensional identity matrix.
GPLVM independently draws a latent coordinate xi from N (0, IQ ) anddraws the dth column of observations y·d = {yid }N
i=1 from the distribution below:p(y ·d |X) = N (y ·d |0, K + β −1 I N ),(1)

where K is the Gram matrix, each component is given by the kernel function
k(x, x′ ) evaluated at two coordinates x and x′ , and β > 0 is the precision.
Latent space sometimes has clusters; however, GPLVM fails to capture themas it assumes a unimodal Gaussian distribution as prior to X. iWMM assumes
the infinite Gaussian mixture model (∞-GMM) [4] defined under Dirichlet process theory [27] as prior to X to model latent clusters. Each coordinate xi isdrawn from the following distribution:p(xi ) =∞X

πk N (xi |mk , R−1
k ),

(2)

k=1

where πk , mk , and Rk represent the mixing weight, mean, and precision matrixof the kth Gaussian distribution, respectively. For simplicity, Rk is assumed to
be a diagonal matrix in this study.

3

Proposed Methods

We utilize neural network Gaussian processes (NNGP) [2], which are equivalent
to a neural network whose weights and biases are marginalized, to eliminate
the necessity to optimize neural weights, biases and the width of hidden layer
while implicitly utilizing a neural network. Subsequently, we explain NNGP and
introduce a latent variable model that combines iWMM and NNGP. Finally, we
propose the NPDV.
3.1

Neural Network Gaussian Processes

In an L-layer fully connected neural network, let ϕ(·) and Nℓ denote the activation and width of the ℓth layer, respectively. It is assumed that the weights andbiases of the ℓth layer, Wijℓ and bℓi for i = 1, 2 · · · , Nℓ and j = 1, 2, · · · , Nℓ−1 , areindependently drawn from N (0, σw /Nℓ ) and N (0, σb /Nℓ ), respectively. A preactivation of the ℓth layer, aℓi (xn ) , is then computed as a linear combination of
Nℓ−1
the post-activations of the (ℓ − 1)th layer {ϕ(aℓ−1(xn ))}j=1
:
i
Nℓ−1

aℓi (xn ) = bℓi +X
j=1

Wijℓ ϕ(aℓ−1(xn )) .
j(3)

Nonparametric Bayesian Deep Visualization

5

tmp
Fig. 3: ARD weights estimated by
NPDV on simulation data. The
three dimensions in red are selected
as the necessary dimensions.

Fig. 4: Generative process of NNiWMM. X is drawn from the ∞-GMM,
and Y is then drawn from Gaussian processes with the ARD-NNGP kernel.

aℓi (xn ) is distributed with a Gaussian distribution from the central limit theorem
when Nℓ → ∞ because the summation in (3) is the sum of i.i.d. random variables.
Because this result holds for any n, the output of the ℓth layer {aℓi (xn )}N
n=1 isjointly distributed with the Gaussian process with the Gram matrix K ℓ . Eachcomponent of K ℓ , k ℓ (x, x′ ), is computed as follows:2k ℓ (x, x′ ) = σb2 + σw
Eaℓ−1 ∼GP(0,K ℓ−1 ) [ϕ(aℓ−1(x))ϕ(aℓ−1(x′ ))] .
i
i
i

(4)

The NNGP kernel is computed by iterating the recursion (4) L times, where
Gaussian processes with this kernel are equivalent to an L-layer ∞-width neuralnetwork whose weights and biases are marginalized. The first step of the original
2 T ′NNGP kernel is the inner product of the input: k 0 (x, x′ ) = σb2 + Q−1 σw
x x.
In contrast, we introduce automatic relevance determination (ARD) [3] weights0′
γ = {γq }Q
q=1 into k (x, x ) to estimate the importance of each dimension:2
k 0 (x, x′ ) = σb2 + Q−1 σw

QX

γq xq x′q .

(5)

q=1

An ARD weight γq increases if the qth dimension is highly related to the observations, and becomes zero if it is totally irrelevant. Fig. 3 shows the ARD weights
estimated by NPDV and the selected dimensions when conducting the simulation study . In the simulation, 20 dimensional latent coordinates were estimatedfrom a 100-dimensional data whose intrinsic structure is in three-dimensional
space. We can observe from Fig. 3 that ARD correctly determines the dimensionality of the intrinsic space. Hereafter, the NNGP kernel with ARD weights
is referred to as the ARD-NNGP kernel and is denoted as k L (xi , xj ). The Grammatrix computed using k L (x, x′ ) is denoted as K L .
K L is generally intractable due to the nonlinearity of X; however, it can beanalytically computed when the activation is an identity map or a nonlinear map
belonging to the polynomial rectified nonlinear function family [28], including

6

H. Ishizuka and D. Mochihashi

ReLU. We use kernel functions compatible with an identity map or ReLU tocompute K L . Appendix A presents the concrete forms 4 .
3.2

NN-iWMM

NN-iWMM is formulated by substituting the ARD-NNGP kernel into iWMM
and infers the latent coordinates, their cluster assignments and the number oflatent clusters. Fig. 4 shows its generative process consisting of generating latentcoordinates X from ∞-GMM and mapping X to the observation space using
Gaussian processes with the ARD-NNGP kernel.
In ∞-GMM, the mixing weights {πk }∞k=1 are drawn from the stick-breaking
process GEM(α) [29], which defines a distribution equivalent to the DirichletQk−1
process. Specifically, πk is computed as πk = ψk j=1 (1 − ψj ), where each ψj
is drawn from the beta distribution Beta(1, α). Subsequently, the mean mk anddiagonal components of the precision matrix {rkq }Q
q=1 of the kth Gaussian distribution are drawn from N (0, I) and the gamma distribution Gam(1, 1), respectively.
For i = 1, 2, · · · , N , the cluster assignment zi is drawn from the categorical
distribution Cat({πk }∞k=1 ) and xi is generated from the zi th Gaussian distribution N (mzi , Rzi ). Then, the Gram matrix K L is computed from X using (4)and (5). y·d is drawn from N (0, K L ) for d = 1, 2, · · · , D.
Due to the absence of weights and biases, it is unnecessary for NN-iWMM
to optimize the weights and layer widths while leveraging a neural network.
Additionally, unlike neural clustering models, using ∞-GMM and ARD weights
allows to determine the number of latent dimensions and clusters, which are
critical to model the latent space appropriately, without tuning.
3.3

Nonparametric Bayesian Deep Visualization

NPDV jointly estimates X and the visual coordinates V by integrating the two
reduction steps, as shown in Fig. 2. We use weighted latent coordinates Xγ T as
the input to prioritize the necessary dimensions of X in terms of ARD weights
when estimating V . Denoting RDR (Xγ T , V ) and L(Y , X) as the visualization
loss to estimate V and the loss of NN-iWMM to estimate X, respectively, we introduce how to integrate these two losses into a Bayesian model using regularized
Bayesian inference (regBayes) [30].
The regBayes framework enables the design of Bayesian models while considering the appropriate constraints on its posterior. This framework is built
on a variational formulation of the Bayesian posterior p(θ|Y ) ∝ p(Y |θ)p(θ),
where Y and θ are the observations and parameters, respectively. p(θ|Y ) canbe viewed as a solution to the following variational optimization problem [31]:Z min KL[q(θ)∥p(θ)] − q(θ) log p(Y |θ)dθ
q(θ)
(6)


s.t.
q(θ) ∈ P ,4

All appendices are provided in the Supplemental Materials.

Nonparametric Bayesian Deep Visualization

7

where P is a set of probability distributions. In the regBayes framework, we
consider the following optimization problem with constraints regarding the expectation of the regularizer Eq(θ) [R(θ, Y )]:Z min KL[q(θ)∥p(θ)] − q(θ) log p(Y |θ)dθ
q(θ)
(7)


s.t.
Eq(θ) [R(θ, Y )] ≤ 0, q(θ)∈ P .
The optimal solution of (7) is obtained byq ∗ (θ) ∝ p(Y |θ)p(θ) exp(−λR(θ, Y )) ,(8)

where λ is the Lagrange multiplier. From (8), the optimal posterior of θ that satisfies the constraints is obtained by its right-hand term. For NPDV, p(Y |θ)p(θ)
corresponds to the likelihood of NN-iWMM L(Y , X):∞
L(Y , X) = p(Y |X)p(X|z)p(z|{mk , r k , ψk }∞k=1 ))p({mk , r k , ψk }k=1 ) .

(9)

R(Y, θ) is given by RDR (Xγ T , V ). Therefore, the optimal posterior of the
NPDV is obtained byq ∗ (X, V ) ∝ L(Y , X) × exp(−λRDR (Xγ T , V ) .

(10)

Hence, the joint optimization of L(Y , X) and RDR (Xγ T , V ) is the posterior
inference of q ∗ (X, V ). NPDV makes the posterior of X suitable for visualization
because RDR (Xγ T , V ) in (10) serves as a regularizer to infer the posterior of
NN-iWMM. λ is a hyperparameter that balances L(Y , X) and RDR (Xγ T , V ),and NPDV degenerates to NN-iWMM when λ = 0. V is treated as a deterministic parameter, as are the NNGP parameters {σw .σb , γ}.
The any dimensionality-reduction method can be used for RDR . In this paper
we combine matrix factorization (MF) and t-SNE, which are widely used inmany different domains, with NPDV. We call these methods NPDV(MF) andNPDV(t-SNE), respectively. Notably, λ = N D practically works well for bothmethods.
NPDV(MF): MF approximates a matrix using a product of two low-rankmatrices, and is one of the widely used linear dimensionality-reduction methods.
For NPDV(MF), the weighted coordinates Xγ T are factorized as the rank S
matrix W and visual coordinates
V : Xγ T ≈ W V T . Denoting ∥ · ∥2F as theP
22
Frobenius norm, ∥A∥F ≡ i,j aij , A ∈ RI×J , and the associated visualization
loss RMF (Xγ T , V ) can be computed as
RMF (Xγ T , V ) = ∥Xγ T − W V T ∥2F .

(11)

NPDV(t-SNE): NPDV(MF) may fail to capture the nonlinear pattern as it
linearly reduces the dimensionality of Xγ T . We then introduce NPDV(t-SNE),
which uses the t-SNE loss for RDR (Xγ T , V ) as a nonlinear counterpart. For
NPDV(t-SNE), RDR (Xγ T , V ) is computed as the divergence of similarities in8

H. Ishizuka and D. Mochihashi

the latent and visual spaces. Denoting ⊙ as the element-wise Hadamard product,
the similarity of weighted latent coordinates, xi ⊙ γ and xj ⊙ γ, is computedfrom the conditional probability based on the Gaussian kernel:XpX
exp(−∥γ ⊙ xi − γ ⊙ xj ∥2 /2τi2 )
i|j + pj|i
X
,
p
=
,
ij
2
2
2N
ℓ̸=i exp(−∥γ ⊙ xi − γ ⊙ xℓ ∥ /2τi )pX
j|i = P

(12)

where τi2 is the variance of the Gaussian distribution and is computed fromthe neighbors of γ ⊙ xi using a binary search with perplexity ρ, which is a
hyperparameter that controls the number of neighbors.
The similarity of two visual coordinates, vi and vj , is evaluated using Student’s t-distribution kernel as follows:(1 + ∥vj − vi ∥2 )−1.
pVij ≡ P P
2 −1kℓ̸=k (1 + ∥vk − vℓ ∥ )(13)

VRt-SNE (Xγ T , V ) is computed as the KL divergence between {pXij }i,j and {pij }i,j :Rt-SNE (Xγ T , V ) ≡X
i,j,i̸=j4pX
ij log

pYij
.
pVij

(14)

Bayesian Training

We employ variational inference to train the NPDV. NN-iWMM is a special caseof NPDV when λ = 0; hence, we focus on the training algorithm for NPDV. The
parameters of NPDV are estimated by maximizing the evidence lower bound(ELBO) L of NPDV. L is derived from Jensen’s inequality and the variational
distribution Qthat is used to approximate the true posterior [32]:L(Y , X) × exp(−λRDR (Xγ T , V )L = EQ log.
(15)
Q
Following [33], the ∞-GMM is approximated by a finite Gaussian mixture model
whose maximum number of mixtures is K. Additionally, we impose the following
mean-field assumption on Q.
Q=

N
Y
i=1

q(xi )q(zi )KYq(ψk )q(mk )q(r k ) ,

(16)

k=1

where q(xi ) = N (µi , Si ). The variational distributions for the mixture modelq(ψk ), q(mk ), q(r k ), and q(zi ) have the same form as in [33]. Using the mean-fieldassumption on Q and conditional independence of NN-iWMM, L is decomposed
into the four terms:
L = Eq(X) [log p(Y |X)] − Eq(X) [log q(X)].
p(X, z, {mk , r k , πk }K
k=1 )+ Eq(X,z,m,r,ϕ) log− λEq(X) [RDR (Xγ T , V )] (17)q(z, {mk , r k , πk }K
k=1 )PN= L1 + i=1 H(q(xi )) + L2 − λRDR (Xγ T , V ), .

Nonparametric Bayesian Deep Visualization

9

Algorithm 1 Variational inference algorithm for NPDV
Input observations Y and the number of layers L
1. Pre-training
2
2
Initialize Π 0 = [{µi , S i }N
i=1 , ζ, σb , σw , γ, β]for i=1,2,... do
Update Π 0 with a gradient-based method
end forInitialize V
2. Training NPDV.
Initialize Π 1 = {z, {mk , r k , πk }K
k=1 }for i=1,2,... dof with (18)Generate X
fApproximate L1 , L2 , and RDR with X
Update Π 1 with the EM algorithm in Appendix B22
Update Π 2 = {µi , S i }N
i=1 , ζ, V , σb , σw , γ, β} with a gradient-based method
end forUnlike Gaussian entropy H(q(xi )), L1 , L2 , and RDR (Xγ T , V ) cannot be evaluated analytically. Therefore, we adapt the reparameterization trick [34] to approximate these quantities. Using this trick, Monte Carlo samples of xi , xei , isgenerated by affine-transforming the standard Gaussian noise ϵ with µi and Si :ei = µi + S i ϵi ; ϵi ∼ N (0, I Q ) for i = 1, 2, · · · , N .
x(18)

f = {eHereafter, we outline the evaluation of L1 , L2 , and RDR using X
xi }N
i=1 .
3
The naive approximation of L1 requires O(N ) complexity because it inverts K L ∈ RN ×N and is difficult to train over a large dataset. We exploit
the inducing-point approach [35] to reduce the complexity. Using this approach,f is approximated by O(M 3 ), (M ≪ N ) complexity based on M
L1 given X
pseudo-inputs ζ ∈ RM ×Q in the latent space, and the corresponding Gaussian
process outputs ud ∈ RM . The cost of the variational mixture model L2 , given
f is of the same form as the ELBO for ∞-GMM in [33]. Moreover, the paramX,
eters of the variational mixture model can be updated using the expectationf
maximization (EM) algorithm. RDR (Xγ T , V ) is computed by substituting X
and ARD weights with γ. Appendix B provides the details of evaluating L1 andL2 and the update formulae of the parameters of the variational mixtures.
Algorithm 1 summarizes the NPDV training algorithm. First, we pretrain
the NN-iWMM that assumes N (0, IQ ) as the prior of X to initialize Π0 =22[{µi , S i }N
i=1 , ζ, σw , σb , γ, β]. Subsequently, we initialize V . After that, We generf using (18) and approximate L1 , L2 , and RDR to train
ate Monte Carlo samples X
the NPDV. Then, we update the parameters of the variational mixture modelN22
Π1 = [z, {mk , r k , ψk }Kk=1 ] and the others Π2 = [{µi , S i }i=1 , ζ, V , σb , σw , γ, β]using the EM algorithm and a gradient-based method, respectively.
Less parametricity: Due to the absence of weights and biases, NPDV hassignificantly less parameters than autoencoder variants. A symmetric neural au-

10

H. Ishizuka and D. Mochihashi

Fig. 5: 3D Visualization of a 100-dimensional data that are generated by transforming the mammoth data in (a) using a neural network. (b)–(d) are obtainedby proposed methods and colored by the estimated clusters. (e)–(j) are obtainedby existing methods,
PL
toencoder has 2 ℓ=1 (Nℓ + Nℓ Nℓ−1 ) weights and biases, where Nℓ and N0 represent the width of ℓth hidden layer and dimensionality of the observations,
respectively. It is not uncommon for them to have more than 107 parameters as
Nℓ often exceeds thousands. Conversely, because NNGP has no need to estimate
weights and biases, NPDV has much less parameters than neural models.
Hyperparameter settings: NPDV has several hyperparameters; however,
not all of them are considered in practice. For the maximum latent dimensions
and cluster, Q and K, if they are set to sufficiently large values, the necessary
latent dimensions and number of clusters are estimated by the ARD mechanism
and ∞-GMM. For the perplexity of NPDV(t-SNE), ρ, the number of inducing
points M , learning rate η, and balance term λ, NPDV(t-SNE) achieves higher
accuracy than the existing methods with M = 100, ρ = 30, η = 0.01, and λ = N D.
Therefore, we only focus on a single hyperparameter, the layer depth of the
NN-iWMM L.

5

Simulation Study

We present the qualitative properties of the NPDV through a simulation study.
As mentioned in Section 1, the visualization accuracy may degrade when observations are distributed on a lower dimensional manifold. To imitate such situation,
the 3D mammoth data in Fig. 5(a) is embedded into a 100-dimensional space by
a neural network. Then, we visualize this 100-dimensional data in 3D space using
several methods. Besides NN-iWMM, NPDV(MF) and NPDV(t-SNE), we apply
six existing methods to the data as baselines: Autoencoder5 , t-SNE, UMAP [10],Trimap [14], NCvis [9], and PaCMAP [15]. Appendix C provides the details of
the experimental settings.
5

The network architecture is the same as that of the data generation network

Nonparametric Bayesian Deep Visualization

11

Table 1: Summary of the datasets used in Section 6. C and D are the number oflabels and dimensionality of observations, respectively. For 20 news, 20 labels
are converted to 6 meta-labels according to the dataset guideline.
type
DatasetMNISTImages
Fashion-MNISTImagesDocuments
20 news

C
10
10
6

D
784
784
1,000

Table 2: Average k-nearest neighbors classification accuracy with five different
random seeds. The highest scores are in bold font.
MNIST
Fashion-MNIST
20 news
Methodk=10 k=20 k=30 k=10 k=20 k=30 k=10 k=20 k=30
t-SNE
0.930 0.920 0.915 0.819 0.794 0.783 0.726 0.700 0.690
UMAP
0.921 0.916 0.913 0.777 0.763 0.757 0.729 0.705 0.695
Trimap0.902 0.897 0.891 0.774 0.760 0.757 0.740 0.720 0.713
0.891 0.886 0.884 0.776 0.764 0.759 0.405 0.358 0.338
NCVis
PaCMAP
0.902 0.896 0.894 0.778 0.766 0.759 0.741 0.724 0.717
VSB-DVM
0.931 0.920 0.915 0.837 0.819 0.806 0.778 0.757 0.749
NN-iWMM0.893 0.884 0.881 0.765 0.748 0.741 0.725 0.706 0.698
0.529 0.495 0.484 0.650 0.629 0.619 0.636 0.618 0.610
NPDV(MF)
NPDV(t-SNE) 0.928 0.917 0.911 0.834 0.820 0.808 0.786 0.761 0.750

Fig. 5 shows the resulting plot of each method. In contrast to conventional
visualization methods in Fig. 5 (f)–(j), since NN-iWMM based methods infers
clusters in addition to latent coordinates, we colored the associated plots by
estimated clusters. Especially, Fig.5 (b) shows that NPDV(MF) recovers theoriginal mammoth shape more accurately than the existing methods in Fig.5(e)—(j) and other NN-iWMM based methods in Fig.5 (c) and (d). This means
NPDV(MF) accurately recovered the intrinsic manifold embedded in a highdimensional space in this simulation. Furthermore, it enables to find the specificparts of the mammoth body, such as paws and horns, as clusters by coloringwith the cluster assignments.
The plot of NPDV(t-SNE) in Fig.5 (c) is blurred because the t-SNE occasionally fails to capture the global structure. However, as shown in the next section,
NPDV(t-SNE) achieves superior performance to NPDV(MF) on real-world data.

6

Experiments on Real-World Data

We demonstrate several advantages of NPDV(t-SNE) through real-world dataexperiments on three datasets. MNIST contains hand-written digit images,
where each image is labeled one of 0–9. Fashion-MNIST contains images ofclothing, where each image is labeled with one of 10 categories, such as Tshirts or shoes. 20 news corpus records English articles, where each article isclassified into one of 20 labels. Table 1 summarises these datasets. We randomly

12

H. Ishizuka and D. Mochihashi

Fig. 6: Visualization of Fashion-MNIST. (a), (b) and (c) are colored by latentclusters that can be discovered by each method.

Fig. 7: Visualization of 20 news. (a), (b) and (c) are colored by latent clusters
that can be discovered by each method.
extracted 5,000 samples from each dataset for evaluation. For MNIST and
Fashion-MNIST, all images are scaled within [0, 1]. For 20 news, the original
20 labels were converted into six meta-labels, e.g., rec and sci for recreation andscience, respectively, according to the dataset guideline 6 , and documents weretransformed into 1,000-dimensional tf.idf vectors after removing stopwords andperforming lemmatization.
We used the k-nearest neighbor classification accuracy for k = [10, 20, 30] asa metric. This metric increases when coordinates with the same label are close
to one another and measures how accurately they can capture label differences.
Furthermore, we qualitatively compared the resulting plots. For model comparison, in addition to the methods used in section 5, we built VSB-DVM+t-SNE
that applies t-SNE to latent coordinates estimated by the latest neural clustering model, VSB-DVM [23]. The hyperparameters of VSB-DVM are tuned byminimizing the loss of held-out 1,000 samples using Optuna [36] with 50 trials.
Note that tuning VSB-DVM is time-consuming due to iterative model fitting.
For NPDV(MF) and NPDV(t-SNE), the layer depth, L, the maximum numberof latent dimensions and clusters, Q and K are set to L = 6, Q = 100 and
K = 50, respectively. Appendix D provides The details of experimental settings
and all visualization results.
Quantitative and qualitative comparison Table 2 lists the k-nearest neighborclassification accuracies. VSB-DVM+t-SNE and NPDV(t-SNE) outperform othermethods for Fashion-MNIST and 20 news. Notably, NPDV(t-SNE) shows
comparable accuracies with well-tuned VSB-DVM+t-SNE. Fig. 5 and 6 show the6

http://qwone.com/˜jason/20Newsgroups/

Nonparametric Bayesian Deep Visualization

13

Table 3: Number of parameters of optimized models.
MethodVSB-DVM+t-SNE
NPDV(t-SNE)

MNIST
133.6 × 106
1.04 × 106

Fashion-MNIST
29.1 × 106
1.04 × 106

20 news
2.81 × 106
1.04 × 106

Table 4: Elapsed time for tuning. h and m is hours and minutes, respectively.
MethodVSB-DVM+t-SNE
NPDV(t-SNE)

MNIST
5 days 13h 47m
8h 32mFashion-MNIST
20 news
4 days 23h 23m 6 days 0h 13m
8h 6m
8h 19mresulting plots on Fashion-MNIST and 20 news obtained by the five methods.
For NPDV(t-SNE), NN-iWMM and VSB-DVM+t-SNE, the points are coloredby the estimated clusters. The coloring helps to understand cluster structures
more easily than t-SNE and PaCMAP, as in the simulation study. For 20 news,
NPDV(t-SNE) shows better cluster separation compared to VSB-DVM+t-SNE.
We guess the reason why cluster structures are taken over to visual coordinates
due to the joint training of NN-iWMM and t-SNE.
Computational cost comparison with neural clustering model NPDV(tSNE) incurs significantly less computational cost compared to VSB-DVM+tSNE. VSB-DVM needs to estimate numerous weights. Consequently, it has
2.8–133 times more parameters than NPDV(t-SNE), as shown in Table 3. Additionally, VSB-DVM must tune several hyperparameters, which increases thecomputational time. Specifically, it took multiple days with 50 trials, as shownin Table 4. Conversely, we only focus on layer depth L to train NPDV(t-SNE).
NPDV(t-SNE) finishes computation in a considerably shorter time than VSBDVM even if we try all candidates L = {4, 5, 6, 7}.
Latent cluster discovery We investigated the characteristics of the clusters estimated using the NPDV(t-SNE). In addition to the 20 news, we used
the Brown corpus 7 . For the Brown corpus, we randomly extracted 5,000
sentences and converted them into sentence vectors using SIF weighting [37].
For 20 news, we investigated therelationships between the four lower
labels in rec and estimated clusters.
The articles in rec belong to one of
23 clusters. Fig. 8 is a cross table ofthese labels and clusters. Evidently,
the same labeled articles belong to
specific clusters, while these labels
were not used during training. The
link between labels and clusters in- Fig. 8: Cross table of lower labels in recdicates that a common topic exists (x-axis) and the clusters (y-axis).
7

http://korpus.uib.no/icame/manuals/BROWN/INDEX.HTML14

H. Ishizuka and D. Mochihashi

Fig. 9: Visualization of Brown corpus using NPDV(t-SNE). The points are colored by the estimated clusters. The right is an entire plot, and the left shows
three clusters and their sample sentences.
within a cluster because each label represents a single topic. Fig. 9 shows a scatter plot of the Brown corpus. We discovered that some clusters shared topics,
such as religion, social thought, and science. Furthermore, similar themes (religion and social thought) were placed close, whereas different themes (science)
were distant, based on the first two themes. Therefore, the distance between
clusters reflects the similarity of themes.
From these visualizations, we found that a common topic exists within acluster and the distance between clusters reflects the similarity of topics. As
topics can be considered as intrinsic clusters in a dataset, NPDV(t-SNE) can
help to reveal clusters and grasp their similarities.

7

Conclusion

We proposed a nonparametric Bayesian latent variable model, NN-iWMM, andan associated visualization method, NPDV. NN-iWMM determines the layer
widths, the dimensionality of latent space, and the number of clusters that are
critical to model the latent space without tuning, while leveraging the power ofneural networks implicitly. NPDV estimates the optimal latent coordinates tolearn visual coordinates by integrating NN-iWMM and a visualization method.
Additionally, we introduced NPDV(MF) and NPDV(t-SNE). Both methods enable to visualize the internal structure of dataset by utilizing the estimatedclusters. Simulation studies demonstrated that NPDV(MF) infers the intrinsic latent manifold better than the existing methods. Real data experiments
demonstrated that NPDV(t-SNE) outperforms conventional methods and shows
comparable accuracy with a well-tuned neural clustering model. Furthermore,
it shows two preferable properties in unsupervised settings: (1) NPDV(t-SNE)
takes considerably less training time than the neural clustering model and (2) it
has the ability to revealing plausible clusters without label information..
In this paper, we limit ourselves to study the properties of NPDV in the caseof using matrix factorization or t-SNE. By considering combination with othermethods, we expect to improve the accuracy and computing efficiency of NPDV.

Nonparametric Bayesian Deep Visualization

15

References
1. van der Maaten and Hinton, G.: Visualizing data using t-SNE, Journal of Machine
Learning Research 9(1), 2579–2605, (2008)
2. Lee, J. H., Bahri, Y., Novak, R., Schoenholz, S. , Pennington, J. and Sohl-Dickstein,J.: Deep neural networks as gaussian processes, International Conference on Learning Representation 2018 48 478–487 (2018)
3. Mackay, D. J. C: Bayesian Non-Linear Modeling for the Prediction Competition,
ASHRAE Transaction100(2) 1053–1062 (1994)
4. Rassmusen, C. E.: The infinite Gaussian mixture model, Proceedings of the 12th
International Conference on Neural Information Processing Systems, 554–560.
5. Kruscal, J. B.: Multidimensional scaling by optimizing goodness of git to a nonmetrichypothesis, Psychometrika 29 1–27 (1964)
6. Tenebaum, J. B., de Silva, J. C., Langford, A: A global geometric framework for
nonlinear dimensionality reduction, Science 290(5500) 2319–2323 (2000)
7. Hinto, G. E. and Roweis S.:Stochastic Neighbor Embedding, Advances in Neural
Information Processing Systems 15 (2002)
8. Tang, J., Liu, J., Zhang, M. and Mei, Q.: Visualizing large-scale and highdimensional
data. the 25th International Conference on the World Wide Web 287-–297 (2016)
9. Aleksandr, A. and Maxim, P.: NCVis: Noise Contrastive Approach for Scalable
Visualization, Proceedings of The Web Conference 2020 2941–2947 (2020)
10. McInnes, L., Healy, J., Saul, N. and Großberge, L:UMAP : Uniform Manifold
Approximation and Projection, The Journal of Open Source Software 3(29), 2579–
2605 (2018)
11. Hadsell, R., Chopra, S. and LeCun, Y: Dimensionality reduction by learning aninvariant mapping. In Proceedings of the IEEE Computer Society Conference on
Computer Vision and Pattern Recognition 2 1735-–1742 (2006)
12. van der Maaten and Weinberger, K.: Stochastic triplet embedding. In Proceedingsof the IEEE International Workshop on Machine Learning for Signal Processing,
1-–6 (2012).
13. Wilber, M.J., Kwak, I. S., Kriegman D. J., and Belongie, S.] Learning concept
embeddings with combined human-machine expertise. In Proceedings of the IEEE
International Conference on Computer Vision2, pages 981—989. (2015).
14. Ehsan, A. and Manfred K. W: TriMap: Large-scale Dimensionality Reduction Using Triplets, arXiv preprint arXiv:1910.00204 (2019)
15. Wang, Y., Huang, H. M., Rudin, C. and Shaposhnik, Y.: Understanding How
Dimension Reduction Tools Work: An Empirical Approach to Deciphering t-SNE,
UMAP, TriMap, and PaCMAP for Data Visualization, Journal of Machine Learning
Research 2021(201) 1–73 (2021)
16. Wallach, I. and Liliean, R. : The Protein-Small-Molecule Database, A NonRedundant Structural Resource for the Analysis of Protein-Ligand Binding, Bioinformatics25(5) 615-–620 (2010)
17. Hamel, P. and Eck, D.: Learning Features from Music Audio with Deep Belief
Networks, Proceedings of the International Society for Music Information Retrieval
Conference: 339-–344 (2010)
18. Geng, X.; Zhan, De-Chuan. and Zhou, Zhi-Hua: Supervised nonlinear dimensionality reduction for visualization and classification, IEEE Transactions on Systems,
Man, and Cybernetics 35(6) 1098–1107 (2005)
19. Venna, A.; Peltonen, J.; Nybo, K.; Aidos, H. and Kaski, S.:Information Retrieval
Perspective to Nonlinear Dimensionality Reduction for Data Visualization, Jounalof Machine Leaninrg Research 11(13) 451–490 (2010)

16

H. Ishizuka and D. Mochihashi

20. Zheng, J.; Hangke Zhang, H.; Cattani, C. and Wang, W.: Dimensionality Reduction by Supervised Neighbor Embedding Using Laplacian Search. Biomedical Signal
Processing and Modeling Complexity of Living Systems 2014 (2014)
21. Xie, J., d Girshick, R. and Farhadi, A.: Unsupervised deep embedding for clustering analysis, Proceedings of the 33rd International Conference on International
Conference on Machine Learning 48 478–487 (2016)
22. Fard, M. N. and Thonet, T. and Gaussier, E.: Deep k-means: Jointly clusteringwith k-means and learning representations, ArXiv:1806.10069, (2018)23. Yang, X., Yan, Y., Huang, K. and Zhang, R.: VSB-DVM: An end-to-end Bayesiannonparametric generalization of deep variational Mixture Model, 2019 IEEE International Conference on Data Mining (2019)24. Iwata, T., Duvenaud, D. and Ghahramani, Z.: Warped mixtures for nonparametric cluster shapes, Proceedings of the Twenty-Ninth Conference on Uncertainty in
Artificial Intelligence 311–320 (2013)
25. Lawrence, N. D.: Probabilistic non-linear principal component analysis with Gaussian process latent variable models, Journal of Machine Learning Research 6 1783–
1816 (2004)
26. Rassmusen, C. E. and Williams C. K. I.: Gaussian Processes for Machine Learning,
The MIT Press, (2006)
27. Ferguson, T. S.: A Bayesian analysis of some nonparametric problems. 1 209–230,(1973)
28. Cho, Y. and Saul, L. K.: Kernel methods for deep learning, Advances in Neural
Information Processing Systems 22, 342–350 (2009)29. J. Sethuraman: Constructive definition of Dirichlet process, Statistical Sinica 4(2)639–650 (1994)
30. Zhu, J., Chen, N. and Xing, E. P.: Bayesian inference with posterior regularization
and applications to infinite latent SVMs, Journal of Machine Learning 15(1) 1799–
1847 (2014)31. Zellner, A.: Optimal information processing and Bayestheorem, American Statistician 42(4) 278–280 (1988)
32. Bishop, C. M.: Pattern recognition and machine learning, Springer (2006)
33. Blei, D. and Jordan, M.: Variational inference for Dirichlet process mixtures, Journal of Bayesian Analysis 1(1) 121–144 (2006)
34. Kingma, D. P. and Welling, M.: Auto-encoding variational Bayes, Proceedings ofthe 2nd International Conference on Learning Representations (2013)
35. Titsias, M. K.; and Lawrence, N. D.: Bayesian Gaussian process latent variable
model, Proceedings of the 13th International Workshop on Artificial Intelligenceand Statistics 9 844-851 (2010)
36. Akiba, T., Sano, S., Yanase, T., Ohta, T. and Koyama, M: Optuna: A nextgeneration hyperparameter optimization framework, Proceedings of the 25th ACM
SIGKDD International Conference on Knowledge Discovery Data Mining, KDD’19 2623–2631 (2019)37. Arora, S., Liang, Y. and Ma, T. A simple but tough-to-beat baseline for sentence
embeddings, International Conference on Learning Representations (2017)