Inverse Kernel Decomposition
Authors: Chengrui Li
Anqi Wu
Venue: Accepted by TMLR ()

Abstract: The state-of-the-art dimensionality reduction approaches largely rely on complicated optimization procedures. On the other hand, closed-form approaches requiring merely eigen-decomposition do not have enough sophistication and nonlinearity. In this paper, we propose a novel nonlinear dimensionality reduction method---Inverse Kernel Decomposition (IKD)---based on an eigen-decomposition of the sample covariance matrix of data. The method is inspired by Gaussian process latent variable models (GPLVMs) and has comparable performance with GPLVMs. To deal with very noisy data with weak correlations, we propose two solutions---blockwise and geodesic---to make use of locally correlated data points and provide better and numerically more stable latent estimations. We use synthetic datasets and four real-world datasets to show that IKD is a better dimensionality reduction method than other eigen-decomposition-based methods, and achieves comparable performance against optimization-based methods with faster running speeds. Open-source IKD implementation in Python can be accessed at \url{https://github.com/JerrySoybean/ikd}.

Published in Transactions on Machine Learning Research (03/2024)
Inverse Kernel Decomposition
Chengrui Li
entichengruiigglech..cm
School of Computational Science & Engineering Georgia Institute of Technology
Anqi Wu
ongioou@gatech.comm com
School of Computational Science & Engineering
Georgia Institute of Technology
Reviewed on OpenReview: https: // openreview. net/ formmo id- #40E7toXpa
Abstract
The state-of-the-art dimensionality reduction approaches largely rely on complicated opti- mization procedures. On the other hand, closed-form approaches requiring merely eigen- decomposition do not have enough sophistication and nonlinearity. In this paper, we propose novel nonlinear dimensionality reduction method- Berthod-Ivveree Kernel Decomposition (IKD) based on an cigen-decomposition of the sample covariance matrix of data. The method is inspired by Gaussian process latent variable models (GPLVMs) and has comparable perfor- mance with GPLVMs. To deal with very noisy data with weak correlations, we propose two cauttonssblockwee and geodesic-to make use of locally correlated data points and pro- vide better and numerically more stable latent estimations. We use synthetic datasets and four real-world datasets to show that IKD is better dimensionality reduction method than other rigen-decomposition- based methods, and achieves comparable performance against optimizatiionbasee methods with faster rumning speeds. Open-source IKD implementation in Python can be accessed al btttps://github cwm Som//erryyybbaa/ikke
1 Introduction
Dimensionality reduction techniques have been widely studied in the machine learning field for many years, with massive applications in latent estimation (Wu et al. 2017 2018), noise reduction (Sheybani & Javidi, 2009), cluster analysis (Bakrania et al. 2020), data visualization (Van der Maaten & Hinton 2008a) and so Torth. The most commonly used method is the principled component analysis (PCA), ), a linear dimensionality reduction approach. It favored thanks to the easy use of one-step eigen-decomposition. Its simple linear assumption, however, restricts its exploitation, especially in highly nonlincar sccnarios. On the other hand, nonlinear dimensionality reduction models, such autoencoders (Kramer 1991), variational autoencoders
(VAE) (Kingma & Welling 2013), t-SNE (Van der Maaten & Hinton 2008b UMAP (MeInnes et al. 2020).
2005) can achieve state-of-the- art (SOTA) performance in terms of finding (sub) oottmol low-dimensional latent and rendering satisfactory downstream analyses (e.g., visualization, prediction, classification). However, all these nonlinear models involve intricate optimization which time-consuming, casy to get stuck in bad local optima, and scnsitive to initialization. In this paper, we propose & novel nonlinear -egenn-dccopoiiiionn- reduction approach that finds low-dimensional latent with closed-fom fooe solution but intricate nonlinearity.
and Gaussian process latent variable models (GPLVMs) (Lawrence) 2003
dimensionality
The proposed method is called Inverse Kernel Decomposition (IKD), inspired by GPLVMs. GPLVMs are probabilistic dimensionality reduction generalive models that use Gaussian Processes (GPs) to find a lower dimensional nonlinear embedding of digh-dimensional data. GPLVM and its many variants have been proposed in various domains (Bui Turner, 2015 Wang et al. 2005 2008 Urtasun et al. 2006 Wu et al., 2017) and proven to be powerful nonlinear dimensionality reduction and latent variable models. However, GPLVMs are highly nonlinear and non-convex due to the GP component, resulting in practical difficulties

---

Published in Transactions on Machine Learning Research (03/2024)
during optimization. By deriving the relationship implied by the kernel function in GPLVMs, IKD solves the GPIVM problem through eigen-decomposition. which could give a more stable latent estimation in a shorter time than the traditional pptimization-basse GPLVM solver (Sec. 3.1.1).
In the experiment section, we compare IKD against four eigen-deompoatiinn-eee based and four optimization- based dimensionality reduction methods using synthetic datasets and four real-world datasets, and we can summarize four contributions of IKD:
As dngen-decompositionn- based method, IKD achieves more reasonable latent representations than other iigen-decomposition-- based methods with better classification accuracy in downstream classification tasks, The rumning time of IKD is on par with other digen-decompesttio--ass methods.
IKD is able to provide competitive performance against some SOTA optimization- methods but at a much faster running speed.
IKD promises a stable and unique optimal solution up to an affine transformation. In contrast, optimization- based methods do not guarantce unique optimal solutions and sometimes are not numeri- cally stable due to the highly nonconvex optimization landscapes. Therefore, IKD can sometimes achieve better latent representations and classification performance than optimization-basee methods like GPLVM and VAE.
When the observation dimensionality is large (i.e. observation data is high-dimensional), a lot of methods have significant drawbacks. For example, t-SNE, UMAP, and VAE encounter the curse of dimensionality problems. The large dimensionality not only leads to longer rumning time but also hurts the dimension- ality reduction performance. In contrast, IKD always obtains improved performance with an increasing observation dimensionality, and claims its absolute superiority when the observation dimensionality is very large.
Note that we are not claiming to propose the best dimensionality reduction approach that beats all other SOTAs. We propose an advanced digen-decomposstioon based method that (1) outperforms other eigen- ducompositiio--assed methods in most synthotic and real-world applications with the same scale of running speeds and (2) reaches a comparable level against other oplimizatioon-bssdd methods but with much faster running speed.
1.1 Related works
As GPLVM's dign-dccomposttion-beed solver, IKD starts from the data-generating process instead of the observed data. Although IKD makes use of eigen-decomposition and kernel functions, it is different from other eigen-decomposition based methods. For example, GPLVM uses a kernel function to form a nonlinear mapping from the embedded latent space to the data space, which is opposite to the use of kernel as in kernel PCA (Schölkopf et al., 1997). Another similar method is Isomap (Tenenbaum et al. 2000), which is a generalized version of multidimensional scaling (MDS) (Borg & Groenen] 2005). From the form of the backbone IKD algorithm derived in Sec. 2.2 IKD obtains its target similarity matrix to be decomposed from data's pairwise correlations, and Isomap obtains its target similarity matrix from data's pairwise generalized distances. However, the initial goals of IKD and Isomap are different. Isomap hopes to place the dataset into a low-dimensional space while preserving pairwise distances as well-equally scaled as possible, but IKD works as an eigen-decompositionn versioned solver for data generalized from GPVLM. In Sec 2 we introduce the generative model GPLVM first, and then convert the problem to an eigen-decomposition problem step- by-step, and finally solve it.
2 Methodology
2.1 Gaussian process latent variable model
Generative model. Let X E RTN be the observed data where T is the number of observations and N is the observation dimensionality of each data vector. Let Z E RTxM denote its associated latent variables where M is the latent dimensionality. Usually, we assume the latent space is lower-dimensional than the original observational space, leading to M < N. For each dimension of X denoted as X BT, Yn e (1,..., N), GPLVM defines a mapping function that maps the latent to the observation which has Gaussian
2

---

Published in Transactions on Machine Learning Research (03/2024)
process (GP) prior. Therefore, given the finite number of observations, can write
X Lid WOO Vn E NN.
(1)
where K is a T x covariance matrix generated by evaluating the kernel function it GP at all pairs of rows in Z, i.., kjj = k(z,, 2j) where 2, and 2j are the 7th and jth rows of Z.
Problem setting. The goal of GPIVM to estimate the unknown latent variables Z that are used for constructing the covariance matrix K, from the observations X. Note that we only consider noiseless GP for the derivation of IKD, but IKD can deal with noisy observations empirically.
2.2 Inverse kernel decomposition
In this section, we derive a novel nonlinear decomposition method, inverse kernel decomposition (IKD); inspired by GPVLM. Previous work has been solving GPLVM by maximizing the log-likelihood obtain Z from X in onc-step fashion. Now let break this process into two steps: (1) cstimating K from the observalions X, and (2) identifying the latent variables Z from the estimated covariance matrix K. The first step can be solved by estimating K with the unbiased estimator, i.e., sample covariance S = N (X- x1") (x X17)"= N N XXT, where X= N E X ought to be since X.,, are i.i.d. samples from a zero-mean Gaussian (Eq. D.
Therefore, our main focus is to estimate the latent Z given S in the second step. In the following, we focus the discussion of commonly used stationary kernel, the squared exponential (SE) kernel. We will show in Sec. 2.3 that IKD can also work with various stationary kernels.
The SE kernel is defined as k(z,,, zz) = o2 exp (- ||==:= ) where a² is the marginal variance and is the 2 length-scale. Note that a² k(z,,, 2,) ki,s, Vi (T) Let f be the scalar function mapping the scaled squared distance d1 " between latent zz) to the scalar covariance ke,j i.e., k.j (dd,g) = x² exp (- -222 Let us assume we know the true K for now. Since f()) is strictly monotonic, we dii))
can obtain di = 1 (ki.)) 21 ( Writing di,j the matrix form D = (di,) ITTT we have
0 (z2 - 21)² T(22 - 21)
(z) - 22) 22)
(z) 27) )T))(Z) - zT)
1
0
(z2 2T) T(22 (22 - 2T)
D = 12
(2)
(zT - 21) T 21) (zT - x2) ((2T x2)
(K) = [5 (kis)) TxX
where f maps K to D element-wisely We define == = with E = O. Now we have
d.j = (2) 2z) - -2,) z))=p((e [(2 21) - (zj - 21)}= [zz - -2) -1)- -22 -1))
(3)
= = " -ē,) = + +2] - -2ē7'zz.
Since 2 = 0, we have d1,j 2 + E3 - -22t'zz = 2T di Vj E {I,..,, T}. Therefore,
we arrive at an expression of 2 23 as ET² =1( + d1,j diz). Note that dı, d., because of the IT RRTx we could write the matrix
symmetric property. Denote Z [2,, ",, ,2 20 E2
form _ = ( 2,) rxr
0
0
0 [ + di. -
d2.1
Zē =
0
3
(da,1 + d1,2 - d3.2)
1 (da, + d1. d3,T)
= =
(4)
0 (dr,1 + d1,2 - dr.2)
dr
3

---

Published in Transactions on Machine Learning Research (03/2024)
which is Pank-M symmetric positive semi-definite matrix given M < ZZT. Then, Eq. has the unique "reduced" eigen-decompositionn
T. g is the function mapping D to
(AAATTT(VUU...) . VA.UU., (VMU....... (.AAU.........UUUU...
(k) UAU (VA.U.......
(5)
where U.m 10, mm us........
mm
TT RR' is the m column of U E RRTM and A A= didg(AA) 22 Am) dr 0. Note that the unique "reduced" singular value
with A A2
AM AM+1
decomposition of Z is
= - = 2, = IVU + 21, WE (1....... T)
(6)
where 21 represents the reference translation, the length-scale 1 is scaling factor, and V is orthogonal matrix that is responsible for the corresponding rotation and reflection. Since Z and Ú span the same column space such that
||| 1ii-zzj1 22
IIV. - IVUT
k(z;; =0² exp
exp
0² exp
212
2/2
(7)
I Úن - = ki= (Ü Úg:).
= exp
Ú contains all of the low-dimensional information of Z. Therefore, we consider Ú as an estimator of Z. Eq. 5 and Eq. summarize the inverse relationship from a GP kernel covariance matrix to the latent variable.
To date, we are able find the exact estimation Ü given the true GP covariance kernel that is constructed from Z (Eq. 7 In practice, we only have the sample covariance estimator S, and neither rank-M nor positive semi-definite is guaranteed for (S)), Therefore, we try to find its optimal rank-M positive semi-definite approximation, i.e.
minimize XERTXM
(8)
Dax et al. 2014) shows that
Ú= (VA.U......... V
(9)
the optimal solution for any unitarily invariant matrix norm I| where A,,..., Am are the first M largest positive eigenvalues of gS)) and U. U,M are the corresponding eigenvectors. The goodness of the estimated latent via the target loss function Eg. can be quantified by the explained variance ratio
is
X3
We summarize the IKD algorithm in Alg.
Algorithm Inverse kernel decomposition
function IKD(X X E TNNN f)
St Nー X
X17) (x "x1)
serves an estimator of the covariance K
o² < 1 Σ 8ú,
estimate σ² through statistic of the diagonal of S
D +
D serves as an cstimation (Eq. 2
4:
5:
U,A +- eigen-decomposition of g D
Eq. 5
6: Form the optimal latent solution Ú using U and A
Eqq 5
7: return Ú
end function
2.3 IKD with general stationary kernels
Apart from the SE kernel, IKD also works for most commonly used stationary kernels, as long as the kernel function f is invertible (i.e., of strictly monotonic over oc)) and we can find a unique non-negative solution for d (kk.. We summarize the kernels in Tab.
4

---

Published in Transactions on Machine Learning Research (03/2024)
Table 1: Stationary kernels that can be applied IKD.
kernel
f-1
squared exponential
f(d) = 0² exp
(- 11
[-) = -2
rational quadratic
f(d) = o² (1+
za
)
d
f (k) = -a
(
7-exponential
[(d) = oeexpl Pexp(-e (-d2)
(k) = (-Inm no closed-form but solvable with root-finding algorithms
f(d) = o T(v) (Venaa) Ky (V2ivd)
Matérn
For the SE kernel, we can generalize it to the ARD kernel k(z,, 2z) = o² exp Em F. (2 - =jmm)² and the Gaussian kernel k(z,,z) = x² exp (-2 - 2z) T (zz 2,)), where an extra affine transformation L is needed, rather than constant scaling 1.
For the Matérn kernel parameterized by v, (-) ) is the modified Bessel function of the second kind. Although it is complicated to obtain closed-form of f (-) for the Matérn kernel, f-) always exists since f(-) is strictly monotonically decreasing over [0, 00) for all v > 0. Note that for the commonly used v = ppp 1 2. p E N,
it is easy to derive f'(,), e.g., when v = 3 f(d = o² [ v/d
.ad and f'(d) 3do² p- exp 33d
exp
In such cases, higher-order root-finding algorithms (e.g., Newton's method) can be used to solve d = f (k).
The intuition of IKD is to find a non-linear mapping that makes IG, H that are strongly correlated in obsorvational space to be closely located with each other in the latent space. This is perfectly reflected by the shape of these stationary kernels, i.e, monotonically decreasing fast near = and then becoming flat (converge to 0) when d - co. Given these stationary kernels have such a similar shape, the latents estimated by different kernels are also similar (see the last paragraph of the real-world experiment, Sec. 3.2). Since these kernels are all invertible, we can derive the exact relationship between latents (1) and Ú from two kernels fi and fz respectively. Denote D = [(000] - ( 0 TXT
and similar for D(2), then
f1 (Ba)) S = (D!2))
(10)
So, the transformation from ((1) to (22) is the diffeomorphism fz f..
2.4 Error analysis of IKD
IKD performs eigen-decomposition on (S) (Eq. which uses the sample covariance S as an empirical estimator of K. In practice, sample covariance values $i S can be very noisy due to the noise the data and an insufficient observation dimensionality N. There can be non-positive and close-to-zero positive covariance values preventing from calculating dij (str) accurately. Non-positive sij falls out of the input range of f-1 ie., o"). For close-to-zero positive Si,j, the error between the estimation d. = f- and the ground truth dej = f- (kk,)) can be large and sensitive to stj. sketch analysis of the error for the SE kernel is via the Taylor expansion of f-1 at si.i
8hd di (kii)) = -2 $ij + (Ki std)) -2
2 kig - Sad + O((ki.i s1,))))
o²
81.j
(1I)
drs + O(ki, si,)) (sUs) + O si.s)
84.1
8i.
We define the estimation error as |da, dryl -1 For large si,j the error is small; but for small 865
Sigj the error is very sensitive to the covariance error |hi - s.,l. To resolve the issue, there are two solutions:
Blockwise solution. We first throw away bad st,j values by thresholding the sample covariance with a value so, leading to a thresholded covariance matrix $ = (8t. 1 [ sol) $ is not fully connected
5

---

Published in Transactions on Machine Learning Research (03/2024)
graph due to the zero values. We can not directly apply IKD to S for latent estimation. Then, we can use, for example, the Kerbosch algorithm (Bron & Kerbosch, 1973), to find maximal cliques in S. Consequently, each clique is a fully connected subgraph (block) of S, which can be decomposed using IKD. After obtaining the latent for each clique, we merge all of the estimated latent variables according to the shared points between every clique pair. As long as the number of shared points between two cliques is greater than M, we are able to find the unique optimal rigid transformation that aligns every two cliques correctly.
We clarify the details of the clique merging procedure in the blockwise solution here. First, the Bron- Kerbosch algorithm provides us with a set of cliques {Ci, .2...., satisfying U C, = {z) Now, we start with two cliques that share the maximum number of points, arg maxc,o, n,n) For example,
now we have C₁ = {z1, #2, 23, 21} and C2 = {zg, 23, 21, 23} (Fig ! Their shared points are {z2, 23. =1). Since the latent in RM R2 and |C, C, = > M + 1, we can find a unique reflection + rotation + translation via SVD to align {zz, 23. 21} in C2 with {z2, 23, z1} in clique Cj. If there are less than M + 1 shared points between the two cliques, there is more than one way to align them. Similarly, in RRM, the uniqueness of this optimal alignment transformation requires at least M + 1 shared points. Once two cliques are merged, we remove C, Cj from the original set and place the merged cliques C, CJ back. The final latent can be obtained by repeating this merging procedure until the largest clique contains all points.
23
222
2/2
24
255
24
Figure 1: A clique pair example that shares three points (22, 23, 21).
Although the complexity of the Bron-Kerbosch algorithm for finding maximal cliques is O(3²), we can terminate the algorithm as long as the union of the existing cliques is the whole dataset. In other words, we only need up T' maximal cliques, so the clique finding time can be bounded by O(T2). Then solving first M cigen-decomposition algorithm for up to T cliques takes O(T (MT?)). Therefore, the complexity of the entire procedure can be bounded by O(MT8).
Geodesic solution. Since small values si < have significantly bad effects on eigen-decompostion,, we can replace 8s.j, whose value is smaller than sn, with the geodesic covariance 8ij + max (ta.a.......... 8i,t, St..t ... jr where i + t + ................................ + + is the geodesic path from i to j found by the Dijkstra algorithm (Dijkstra et al. 1959.. The complexity of this approach is bounded by the complexity of the Dijkstra algorithm, which is O log(T)). Since the complexity of the geodesic approach is smaller than that of the blockwise approach when T is larger (greater than 1000 in the following experiments), we choose the geodesic instead of the blockwise. A comparison of these two solutions is shown in the experiment section.
2.5 Reference point selection
In Eq. 3 we choose 21 as the reference point to calculate dz,j. But the reference point can be any of the points in {2+]T] 1- Tf we choose 2,, for an arbitrary index r E {1...., T), to be the reference point, then similar to Eqq 4 the th row and the , column of g(D) are Os, and the remaining elements are g(D).,) = (dder dr.j d o O, Vi & r,j # Note that every (D),, includes an element from {dr,)) Thus the quality of [d is vital for latent estimation. Based on the analysis in Eq. 11 we know that in practice we want to choose good reference index r so that {dr 1 are relatively small (i.c., {sr, are large, which means the th data point is highly correlated with the rest of the data points). Note that multidinensional scaling (MDS) (Kruskal & Wish 1978 solves a similar sigen-decompositionn problem, i.e., finding coordinates Z from the distance matrix D. employs a centering idea which is equal to using the average of all latent variables 7 Et= 2, as the reference point. We choose the best reference point instead since we want to reduce the estimation error in Fq. 11 much we can, so that the objective function in Eq. can be
6

---

Published in Transactions on Machine Learning Research (03/2024)
minimized as much possible. Mathematically, we obtain such that ||d, & ||dilbo: Vi f1,2...., T and d: is the 8th row of D = (dı,)rxr, i.c., r = arg min; |d::l>> arg min; max, [dip}}
3 Experiments
In this section, we evaluate IKD with the most commonly used squared exponential as the default kernel on three synthetic datasets, where we know the true latent representations and four real-world datasets.
Baseline methods for comparison:
PCA: Principal component analysis. One of the most widely used lincar dimensionality reduction meth- ods.
KPCA: Kernel PCA. We try different kernels (polynomial, SE, sigmoid, and cosine) and present the best one.
LE (Belkin & Niyogi 2003): Laplacian cigenmaps, which is to do spectral decomposition to the affinity matrix's s graph Laplacian. We use the nearest neighbors algorithm to build the affinity matrix.
Isomap (Tenenbaum et al. 2000): Isometric mapping of multidimensionall scaling (MDS) by incorporat- ing the geodesic distances. We use sklearn's default setting five nearest neighbors to build the distance
matrix.
t-SNE (Van der Maaten & Hinton 2008b): t-distributed stochastic neighbor embedding. We use sk sklearn's default hyperparameter sctting to fit the model.
UMAP (Melnnes et al., 2020): Uniform Manifold Approximation and Projection for dimension reduc-
tion. We use sklearn's 's default hyperparameter setting to fit the model. We use the official UMAP package MMcInnes et al. 2018) and its default hyperparameter sctting to fit the model.
GPLVM (Lawrence) 2003 2005): The traditional optinization-based Gaussian process latent variable model solver. We use the GPLVM module in the GPy package (GPy) since 2012) and its default hyperpa- rameter setting to fit the model. Same IKD, we use the most commonly used squared exponential as the default kernel unless otherwise stated.
VAE (Kingma & Welling 2013): Variational autoencoder.
The first four are eigen-decompositoon-beee based methods; the last four are optimization-base methods.
3.1 Synthetic data
Experimental setup. We first test all methods on three synthetic datasets. All the following experiments are based on 50 independent. repeats (trials). For each trial, we generate the true latent variables from
Z ZT, N (0. (0.(6o- 1==1
Vm E (1,
(12)
TxT)
where M the latent dimensionality, varying across different dalasets. Then, we generate the noiseless data from GP, sinusoidal, and Gaussian bump mapping functions respectively. Afterward, i.i.d. Gaussian noise is added to form the final noisy observations X.
Evaluateion. We evaluate the performance using the R² metric. When computing 1²² values, we first align the estimated latent with the ground truth through an affine transformation (i.e., a linear decoder); then compute R for each latent dimension, and finally take average across all latent dimensions m E (1,2...,My M). The reasons for choosing the affine transformation are: (1) rigid transformation could lead to very negative R²2 values for those non-IKD methods (e.g., PCA), not shown here; and (2) affine transformation is the commonly used one for latent estimation and alignment.
Dataset 1: GP mapping function. We start our experiments with the GP mapping function. In each trial, we generate a 3D latent e g 100003 (i.e., M = 3) according to Eq. 12 and generate X E 10000x according to Eq. [ with o² and 1 Then Gaussian noise addod: Zt,n + Zi,n + Et,n: V(t. n) e (100) {1...., N), where noise εr,, A0 0.053). Note that this generating process is consistent with the generating process of GPIVM. Thus it is well aligned with the model assumptions of IKD, deemed as
7

---

Published in Transactions on Machine Learning Research (03/2024)
(a)
GP
(b)
GP
(c)
GP
1.0
10a
5-
rundie
10a
10 2-
(d)
(e)
2.5
0.0
1.0
107
0.8
0.0
2.5
5.0
7.5
10.0
12.5
15.0
17.5
runtime
0.6
(f)
0.4
2.5
10a²
0.2
-25
(g)
Gaussian bump
(h)
Gaussian bump
0.0
2.5
5.0
7.5
10.0
12.5
15.0
17.5
10a
Gaussian
2.5
rundime
10a
0,0
225
0.4
5
0.2
10:
0
20
100
200
1000
100 200
500 1000
0.0
2.5
5.0
10
12.5 5
N
N
PCA
KPCA
LE
Isomap
1-SNE
UMAP
GPLVM
VAE
IKD
True
Figure 2: R2 values (a,d.g) and rumning time (b,e,h) with respect to N with GP, sinusoidal, and Gaussian bump mapping functions. (c,f,i): Latent recovery visualization of an example trial with N = 100, for the first 20 points with GP, sinusoidal, and Gaussian bump mapping functions.
data-matching example. Fig. 2(a) shows that for N = 10 and N = 20, Isomap is the best; but when N > 50, IKD becomes the best and its 12 converges to as N increases. The latent recovery visualization of an example trial under N 100 for the first 20 points (Fig. [(c) shows that Isomap and IKD match the true latent the best.
Dataset Sinusoidal mapping function. In each trial, we generate 1D latent (i.e., M 1) according to Eq. 12 and generate the noisy observations X 10000XN as a sin( a2n + sp) t&, Vt {1,..., 1000), [PI...., with 9n, u(- R), and noise E, N(0, 0. r).. The result in Fig. d indicates that even though the observed data is not from GP (data- mismatching), IKD is still able to discover the latent structure consistently better than others, except for = 10 and N = 20 where Isomap is the best. When the observation dimensionality N > 50, the R² value of IKD approaches 1, while Isomap, I-SNE, and UMAP all have decreasing performance due to the curse of dimensionality. The latent recovery visualization of an example trial under N 100 for the first 20 points (Fig. (()) shows that Isomap and IKD match the true latent the best.
where S (wn.m)NxM with Wn.m 1),
Dataset 3: Gaussian bump mapping function. In each trial, generate 2D latent (i.e., M 2) according to Eq. 12 and generate X E x as
Ft.n = 20 exp (-11z: c,. + Et...
v(t,n) (1,... 1000} x N1,
(13)
where on 2 is the center of the n Gaussian bump randomly selected from 10,000 grid points uniformly distributed in [-6,00 -]2, with noise ft.x ~ (O, 0.053). This is another data-mismatching example. Fig. 2
shows that in this case, IKD is the best one among all hods for all observation dimensionality N. Fig also shows that only TKD matches the true latent accurately.
The running time of IKD the three synthetic datasets above on par with other digen-dcomposttion-aaaed methods (b,ch) and much less than optimization- based methods. In particular, GPLVM always takes very long time; t-SNE requires more running time datasets whose latent. dimensionality is greater than
8

---

Published in Transactions on Machine Learning Research (03/2024)
squared exponential
rational quadratic
y-exponential
Matérn
1.0
Hard
0.0
A
0.5
-1.0
0.50
Medium
0.25 47
0.00
0
-1
p
1
Easy
R
-0.5
0.5
-0,5
4
Figure 3: Left: ground truth Z of the latent. Right: kernel covariance matrix K created by the corresponding latent for the four diffcrent kernels, with marginal variance a = and length-scale = 0.5; a = in rational quadratic kernel, y = 1 in y-exponential kernel, and v = 1.5 Matérn kernel. The difficulty levels of the three datasets are from hard to easy, from top bottom. For the easy dataset (bottom); points are selected from the grids on the surface and randomly permuted.
2; and VAE is more time-consuminng when the observation dimensionality is large. Note that we only vary the dimensionality N not the number of observations T. When increasing T, the running time of all methods will increase polynomially. For optimization-based methods, stochastie optimization can be employed to scale to large-scale datasets. For fair comparison, extra scaling techniques should be incorporated to deal with large-scale cigen-decompesition, which falls out of the scope of this paper.
In general, IKD performs the best for three mapping functions especially when the observation dimen- sionality N is large. is also very effective in capturing details in addition to recovering the general latent structure correctly. Same as other eigen-decompositionn based methods, IKD takes less time solve the presented problems compared with optimization based methods.
3.1.1 Varying dimensionality, kernels, and latent structures
We test the effectiveness of IKD on the observation data generated from the GP mapping function described in Eq. 1 for different observation dimensionality N (100, 200, 500, 1000, 2000, 5000, 10000), different gen- crating kernels (Tab. D. and three different latent structures (hard, medium, and casy shown Fig. 3 according to their difficulty levels). We only compare IKD with PCA and GPLVM here. PCA is the most commonly used linear method, it serves baseline. GPLVM is the traditional optimization-basee solver, and IKD is our newly proposed cigen-decomposttion-baaee solver for data generated from the GP mapping function.
From Fig. [bb),, we can tell that IKD is always the best for the most commonly used SE kernel. Compared with PCA and GPLVM, IKD is highly effective especially when N is large, where the R² values of IKD are very close
Fig. [(c) shows the latent recovery visualization from one example trial of the medium dataset when N = 1000. The kernel of the generating model is SE. We can tell that TKD matches the ground
9

---

Published in Transactions on Machine Learning Research (03/2024)
(a)
squared exponential
rational quadratic
v-exponential
Matérn
1.0
Hard
0.8
0.6
0.9
Medium
0.8
R²
0.6
0.9
Easy
0.8
0.7
0.6
100
500
2000
10000
100
500
2000
10000
100
500
2000
10000
100
500
2000
10000
N
(b)
PCA
(c)
10a
-1
runtime
10
8
-1
10-1
10-2
20 40
60
80 100
100
200
500
1000
2000
5000
10000
t
Figure 4: (a): R2 values, with respect to N, of PCA, GPIVM, and IKD for different datasets and kernels. (b): Latent recovery visualization of an example trial of the medium dataset with the SE kernel and N 1000. (c): Average running time in seconds (across different kernels, different datasets, and 50 independent trials) of the three methods w.rt. observation dimensionality N. We take averages across different kernels and different datasets because all of them share similar running time results.
truth the best. For the third dimension, particularly, only the estimated latent from IKD reflects the linear increasing trend of the ground truth correctly. In terms of complexity, GPLVM is time-consuming compared with IKD (Fig. Id)). These results indicate that we can use IKD to recover the latent for data generated from the GP mapping function faster and more accurately.
3.1.2 Varying the noise level
To understand the performance of IKD solving GPLVM with different. levels of noise, we use the medium dataset (Fig. 3) and the squared exponential kernel as the generative GPLVM to simulate the noiseless obscrvation by Eq. Then, different leveled noises are added the observation, i.c., Ft,n <- Zt,. + St,ns ). for sd (0,0.1,.... 0.9, 1}. The standard deviation of the Gaussian noise "sd"
where Et,n (00 sd"). represents the noise level.
10

---

Published in Transactions on Machine Learning Research (03/2024)
Fig. 5 shows that IKD is able to have relatively good estimation datasets with low-level noise. The R of IKD gradually decreases as the noise level increases, but it is still better than GPLVM up to the noise GPLVM solver explicitly considers the noise
level of as high as sd = 1. The traditional optimization- term in its algorithm and the performance even slightly improves as the noise increases. Besides, with larger noise, the observation also makes the GPLVM solver run faster.
(a)
R²
(b)
runtime
1.0
10
0.9
0.8
10-1
GPLVM
IKD
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
1.0
noise level
Figure 5: (a): R performances of PCA, GPLVM, and IKD on datasets with different noise levels. (b): The
corresponding running time.
3.1.3 Comparison of the blockwise and the geodesic solutions
To understand the differences between blockwise and geodesic solutions, we apply the blockwise and the geodesic TKD the medium dataset (Fig. 3) respectively and visualize their results. There are two sets of observations, one is from the SE kernel and the other is from the rational quadratic kernel. We only use the SE kernel to solve both of the two observations, constituting model match case and model mismatch case.
b ockwise
geodesic R 0.7172
R² 9878
model match
0.5
0.50
40
100
0.0
0 0
0
100
-0.50
20
40
60
80
100
0.5
1
0.0
1
0.0
-05
0.5
23
25
50
100
12
40
100
R²
R² 0.7006
model mismatch
0.50
50
100
0.50
100
0.25
0.25
0.0
4
5'
0 0
0.25
0.50
25
50
100
0.50
20
40
80
100
0.5
0.5
0.0
4..0
0.5
0
22
100
0
1
13
100
Figure 6: The blockwise (left) and geodesic (right) solutions in model match case (above) and a model mismatch case (bclow). The blue curve the true latent, and the orange curves are estimated latent. In the model mismatch case, although the true tata-generating model is GP with the SE kernel, we still use a simple decoder (i.e., affine transformation) to align the estimated latent with the true latent to compute the R² metric.
In the model match case, the blockwise approach solves the GPLVM, but the geodesic approach does not. Therefore, in the model match case and when the number of data points is not large, the Bron-Kerbosch algorithm can be finished quickly and the blockwise solution should converge the true latent up to an affine transformation (the model match case in Fig.
11

---

Published in Transactions on Machine Learning Research (03/2024)
In the model mismatch case, the tique-mmergnn procedure is vulnerable to numerical errors, and small errors can accumulate during this procedure, especially when the number of cliques is large. And if the number of cliques is large but points are scattered, there will be too many cliques to be found and hence the running time of the Bron-Kerbosch algorithm is unaceeptable. Moreover, we do not know whether the true data- generaling model is GP with the SE kernel (and it very likely that it is not GP, just like in PCA we know the true data-generating model is certainly not simple linear model). In such cases, we treat all data points as a whole and compute their empirical correlation through the geodesic shortest path, and hence using geodesic is more intuitive, globally stable, and efficient (the model mismatch case in Fig. 6
3.2 Real-world data
Dataset. We compare IKD against alternatives four real-world datasets:
Single-coll qPCR (PRC) (Guo et al. 2010;; Normalized mcasurements of 48 genes of single cell at 10
different stages. There are 437 data points in lotal, resulting in X E R R2437x88
Hand written digits (digits) (Dua & Graff 2017|: It consists 1797 grayscale images of hand written digits.
Each one is an 8 x 8 image, resulting in X E 17977x64
COII-20 Nenn al.) [1966): It consists 1440 grayscale photos. For each one of the 20 objects in total, 72
photos were taken from different angles. Each one is 128 x 128 image, resulting in X 11100016384
Fashion MNIST (F-MNIST) (Xiao et al. 2017): It consists of 70000 grayscale images of 10 fashion items (clothing, bags, etc). We use subset of it, resulting in X E 30000784
Evaluation. Since there is no true latent. to compare against, we first estimate the latent. a (2, 5, 10]- dimensional latent space and then use the k-nearest neighbor (R-NN) classifier to evaluate the performance of each dimensionality reduction method. Specifically, we apply 5-fold cross-validation k-NN (A (5, 10, 20}) the estimated {2, 5, 10)-dimensional latent to evaluate the performance of each method on each dataset. The k-NN classification results of different methods under different latent dimensionality M, different datasets, and different choices of k are shown in Fig. aa). Performances of different methods on differ- ent datasets evaluated by the silhouetle score are presented in Fig. 7b).
Results. Comparing TKD with other eigmn-decompesitton-bbeee methods (PCA, KPCA, LF, Isomap), we can conclude that IKD is almost always the best all four datasets, except that when M E (3, 5,13 in the digits datasel, Isomap is better than IKD. When comparing IKD with GPLVM, we find the performances of GPLVM on PCR, digits, and F-MNIST datasets are slightly better than IKD while GPLVM takes too much running time. Specifically, IKD is significantly better than GPIVM on the COIL-20 datasct but only slightly worse than GPLVM on the other three datasets. VAE only performs well on the most complicated dataset F-MNIST, and is much worse than TKD in the other three datasets. Although TKT) is worse than the remaining two optimization-based methods (t-SNE and UMAP), the performance of IKD is the best the COIL-20 dataset. The reason is that the observation dimensionality is very high (N = 16384) in the COII.-20 dataset, and IKD is very effective for high-dimensionall data as shown in the synthetic results (Fig. aa))
2D visualization of the four datasets are shown in Fig. 10 and respectively. Qualitatively, we can see that IKD consistently finds more separate clusters compared with all other oigen-deompositiio--bssss mothods and two optimization-based mothods (GPLVM and VAE) across all four datascts. Therefore, even though IKD is an eigen-decompossiioon- based method, its performance is significantly better than other eigen- decomposition methods on all four real-world datasets, sometimes good as the best optimization- method.
In terms of rumning time (Fig. 7 IKD is on par with Isomap, and these bigen-deeompossiiion-bbbsss methods are significantly faster than those four optimization based methods. Note that if the desired latent dimensionality M > 2, the running time of t-SNE barely acceptable. For the high dimensional COIL-20, the running time values of VAE and GPLVM are extrenely high, getting out of the upper limit of the corresponding axes.
12

---

Published in Transactions on Machine Learning Research (03/2024)
(a)
PCR
digits
COIL-20
F-MNIST
1.0
1.0
0.8
0.9
07
in
0.8
0.8
0.6
0.7
0.6
0.5
1.0
1.0
08
accuracy
0.9
1=10
0.8
0.8
0.7
0.6
0.6
0.7
0.6
1.0
08
0.9
k=20
0.8
0.8
07
0.7
0.6
0.6
2
3
10
5
10
2
3
5
10
(b)
silhouette score
0.4
0.50
0.50
0.1
0.2
0.00
0.0
(c)
102
10a
10
runtime
10-1
10a
10
2
M
5
10
2
10
5
3
10
KPCA
Isomap
t-SNE
Figure 7: (a) k-NN 5-fold cross-validation, (b) silhouette score, and (c) running time, on different methods, different latent dimensionality M, different datasets, and different choices of k.
PCA
KPCA
LE
Isomap
t-SNE
UMAP
GPLVM
VAE
64
2
64 PE
32 ICM
64 TE
TE
Figure 8: Visualization of the dimensionality reduction results of different methods on the PCR dataset.
13

---

Published in Transactions on Machine Learning Research (03/2024)
PCA
KPCA
LE
Isomap
(-SNE
UMAP
GPLVM
VAE
IKD
0
1
Figure 9: Visualization of the dimensionality reduction results of different methods on the digits dataset.
PCA
KPCA
LE
Isomap
t-SNE
UMAP
GPLVM
VAE
IKD
10
15
11
16
12
14
Figure 10: Visualization of the dimensionality reduction results of different methods on the COIL-20 dataset.
PCA
KPCA
Isomap
t-SNE
UMAP
GPLVM
VAE
IKD
1
2
3
Figure 11: Visualization of the dimensionality reduction results of different methods on the F-MNIST dataset.
Dimensionality reduction results of IKD with different kernels. As derived in Eq. 10 in the main content, choosing different kernels does not affect the dimensionality reduction results too much in terms of the latent structure similarity since different kernels have similar shapes. Especially in a classification task, although the exact. relationship between the estimated latent data points under different kernels is not the same, points belonging to the same class consistently locate together in latent space (low dimensional space) no malter what kernel used. Fig. shows that the dimensionality reduction results under different kernels are nearly the same (up to reflection, rotation, and translation). Tab, 2 also shows similar 5-fold cross-validation k-NN classification accuracies under different kernels.
Varying number of points T. Although a theoretical analysis of the time complexity of IKD w.r.t. the number of data points has been provided in Sec. (2.4and discussed in See. 3.1 we are still curious about the comparisons of different methods on datasets varying number of data points T. Therefore, apply
14

---

Published in Transactions on Machine Learning Research (03/2024)
squared exponential
rational quadratic
v-exponential
Figure 12: Dimensionality reduction results of IKD with different kernels for the digits dataset with the latent dimensionality M = 2. in the rational quadratic kernel, and in the 7-exponential kernel.
5-NN 5-fold cross-validation classification accuracies under different kernels and latent dimension- alities M E 3,,3,5 10} for the digits dataset.
Table 2:
k
kernels
2
3
5
10
squared exponential
0.875899
0.85085
0.946049
0.944937
rational quadratic
0.841382
0.821323
0.931574
0.935474
7-exponential
0.837478
0.806288
0.930458
0.933807
squared exponential
0.872006
0.844732
0.936592
0.937696
10
rational quadratic
0.857527
0.825235
0.92212
0.943258
7-exponential
0.854737
0.81242
0.919336
0.93992
squared exponential
0.871453
0.843067
0.928804
0.932683
rational quadratic 0.857521
0.822467
0.906541
0.929341
7-cxponential
0.856408
0.817457
0.908767
0.928231
different methods the F-MNIST dataset with the number of data points T {1000, 2000, 3000, 4000). The k-NN accuracies and silhouette scores in Fig. 13 (a) and (b) show that IKD is consistently the best among all eigen-deconnposiiion-oeee runtime, we can see that IKD faster than optimization- based methods. However, comparing the growth rate of IKD with UMAP and VAE, we can see from the growth rates shown in Fig. 13 ())c that, due to (1) the worse time complexity of IKD and other eigen-decome decomposooooiiionnn based methods and (2) the good scalability of optimization-based methods to large-scale datasets, the impact of the number of data points T to IKD and other eigen-decompositionn based methods gets bigger and bigger T increases.
methods, except for the silhouette score at == 4000. In terms of
4 Discussions
In summary, IKD, cigan-dcoompsiition--assss method, consumes short running time but is able to obtain dimensionality reduction results better than other eigen-deeomposition--eesse methods. When facing high-dimensional observation data, IKD can perform significantly better than all other methods in a very short time.
Note that although eigen-decompssition--aaee methods perform relatively worse than optimization- mothods, the benefit of fast running provides good initialization for sophisticated nonlinear optimization problems, mitigating the numerical instability and multi-modal issues commonly observed in methods such GPLVM and VAE.
15

---

Published in Transactions on Machine Learning Research (03/2024)
(a)
k-NN accuracy
(b)
silhouette score
(c)
runtime
10a
0.175
0.75
0.150
10a
0.125
0.70
0.100
0 0.65
0.075
0.050
10-1
0.025
4000
1000
2000
4000
1000
4000
PCA
KPCA
LE
Isomap
6-SNE
UMAP
GPLVM
VAE
IKD
Figure 13: (a) average k-NN 5-fold cross-validation, (b) silhouette score, and (c) running time, by different methods, on the F-MNIST dataset including different numbers of data points T..
16

---

Published in Transactions on Machine Learning Research (03/2024)
References
Mayur R Bakrania, 1 Jonathan Rae, Andrew P Walsh, Daniel Verscharen, and Andy W Smith. Using di- mensionality reduction and clustering techniques to classify space plasma regimes. Frontiers in Astronomy and Space Sciences, pp. 80, 2020.
Mikhail Belkin and Partha Niyogi. Laplacian Eigenmaps for Dimensionality Reduction and Data Rep-
resentation. Neural Compulation, (5(6):1373-139)) June 2003. ISSN 0899-7667. doi: 10.1162/ 0899766033211780117, Conference Name: Neural Computation.
Ingwer Borg and Patrick JF Groenen. Modern multidimennsonall scaling: Theory and applications. Springer Science Business Media, 2005.
Coen Bron and Joep Kerbosch. Algorithm 457: finding all cliques of undirected graph. Communicaations
of the ACM, 16(9):575-577 September 1973. ISSN 0001-0782, 1557-7317. URL https: acm. wwt..oogg////..oooomm 0tt.//////////////..//////2222.0
doi: 10. 00.1145/3623422362233
Thang Bui and Richard E Turner. Stochastic variational inference for Gaussian process latent variable models using back constraints. In Black Bor Learning and Inference NIPS workshop, 2015.
Achiya Dax et al. Low-rank positive approximants of symmetric matrices. Advances in Linear Algebra &
Matria Theory, 4(03):172, 2014.
Edsger W Dijkstra et al. A note on two problems in connexion with graphs. Numerische mathematik, 1(1): 269 271, 1959.
Dheeru Dua and Casey Graff. UCI machine learning repository, 2017. URL http::/arrchive.iccsseee edu/
GPy. GPy: A gaussian process framework in python. http:///tthibb.co//afffffaadddd since 2012.
Guoji Guo, Mikael Huss, Guo Qing Tong, Chaoyang Wang, Li Li Sun, Neil Clarke, and Paul Robson. Resolution of cell fate decisions revealed by single-cell gene expression analysis from zygote blastocyst. Developmental cell, 18(4))675-685, 2010.
Diederik P Kingma and Max Welling. Auto-encoding variational bayes. arXiv preprint arXiv:3312.6114 2013.
Mark A. Kramer. Nonlinear principal component analysis using autoassociative neural net- works. AIChE Journal, 37(2):233 243, 1991.
ISSN 1547-5905.
doi: 10.1002/aic.
690370209. URL https:// http::///annnnnebarrrrrr birrryyrwol
wiley cow/ http/////abbb/..ooo/ 1002/002/aa aic. 690370209 69000000 6903700000
.print:
Joseph B Kruskal and Myron Wish. Multidimenssional scaling. Number 11. Sage, 1978.
Neil Lawrence. Gaussian Process Latent Variable Models for Visualisation of High Dimensional Data.
In Advances in Neural Information Processing Systems, volume 16. MIT Press, 2003. URL https:// proceedings. wwweeedding..nuripp...
wetee-bbssttttttte html
Neil Lawrence. Probabilistic Non-linear Principal Component Analysis with Gaussian Process Latent Vari-
able Models. Journal of machine learning research, 6(11), 2005.
Leland Mclnnes, John Healy, Nathaniel Saul, and Lukas Grossberger. Umap: Uniform manifold approxima- tion and projection. The Journal of Open Source Software, 3629),861, 2018.
Leland McInnes, John Healy, and James Melville. UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction. arXiv:1802.03./ fes, stalf, September 2020. URL http://arxivi www.aaiwiv..ooaaaaoo aba/ 1802. 03426 aiviv: 1802.03426.
Sameer A Nene, Shree K Nayar, Hiroshi Murase, et al. Columbia object image library (coil-100). 1996.
17

---

Published in Transactions on Machine Learning Research (03/2024)
Bernhard Schölkopf, Alexander Smola, and Klbes-Robert Müller. Kernel principal component analysis. In Wulfram Gerstner, Alain Germond, Martin Hasler, and Jean-Daniel Nicoud (eds.), Arlificial Neural Net- works ICANN'97, Lecture Notes in Computer Science, pp. 583 588, Berlin, Heidelberg, 1997. Springer. ISBN 978-3-540-69620-9. doi: 10.1007 BF0000//B000021..
Ehsan Sheybani and Giti Javidi. Dimensionality reduction and noise removal in wircless sensor network datasets. In 2009 Second International Conference on Compuler and Electrical Engineering, volume 2, Pp. 674-67- 677. IEEE, 2009.
Joshua B Tenenbaum, Vin de Silva, and John C Langford. A global geometric framework for nonlinear dimensionality reduction. science, 290(5500):2319 23119-233333 2000.
R. Urtasun, D.J. Fleet, and P. Fua. 3D People Tracking with Gaussian Process Dynamical Models. In 2006 IEEE Computer Society Conference on Computer Vision and Pattern Recognétion (CVPR volume 1, pp. 238 245, June 2006. doi: N0.119//CVPR22066115 ISSN: 1063-6919.
Laurens Van der Maaten and Geoffrey Hinton. Visualizing data using t-sne. Journal of machine learning research, 9(1)), 2008a.
Laurens Van der Maaten and GeofTrey Hinton. Visualizing data using U-SNE. Journal of machine learning research, 9(11), 2008b.
Jack Wang, Aaron Hertzmann, and David J Fleet. Gaussian Process Dynamical Models. In Advances in Neural Information Processing Systems, volume 18. MIT Press, 2005. URL httts::/ppoceeddigg neurips.cc/
wwtpee///e.a.o-ttaatttt.oo .htm
Jack M. Wang, David J. Fleet, and Aaron Hertzmann. Gaussian Process Dynamical Models for Human Motion. IEEE Transactions on Pattern Analysis and Machine Intelligence, 30(2):283 298,, 2008. ISSN 1939-3539. doi: 10.1109/TPAMM. 207...2007..1166. Conference Name: TEEE Transactions Pattern Analysis and Machine Intelligence.
Anqi Wu, Nicholas A. Roy, Stephen Keeley, and Jonathan W Pillow. Gaussian process based nonlinear latent structure discovery in multivariate spike train data. Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017. URL htttps://proceedinggs wwwpeeiinnnnnnnunnnn.eeee nwwrine..ccccccco wttpppprccppppp/pppppppppp Abstract html
Anqi Wu, Stan Pashkovski, Sandeep Dataa, and Jonathan W Pillow. Learning latent manifold of odor representations from neural responses in piriform cortex. Advances Neural Information Processing Systems, 31, 2018.
Han Xiao, Kashif Rasul, and Roland Vollgraf. Fashion-mnist: novel image datasct for benchmarking machine learning algorithms, 2017.
18


Reviews:
[{"review_id": "4Yme3fMyqZ", "reviewer": "Reviewer_jwfe", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n- The method is simple and easy to implement\n- It has interesting connections with spectral methods\n- It works reasonably well in practice\n\nWeaknesses:\n- Connections to prior standard methods\n- Experiments do not support the different claims\n\nIn my opinion, there are significant connections to standard methods that are missing, especially regarding the proposed method and MDS, which do not seem to be sufficiently discussed (only in 2.5). In fact, the proposed method seems to do almost exactly what MDS would do, with a minor modification. For MDS, one observes measurements of 'distances' (in a broad sense) $\\delta_{ij}$ and seeks to find latent positions of a squared distance matrix (EDM) $d_{ij} = |x_i - x_j |^2$ such that $\\delta_{ij} \\approx d_{ij}$. To reconstruct the points from the set of distances, a very classical approach is based on the eigenvector decomposition of the covariance matrix (see e.g., [Section B, 1] or [3]). In fact, the approach proposed here is exactly the same approach but with a special choice of $\\delta_{ij}$ that comes from a similarity kernel (mapped to a notion of ''distance'' with the $f^{-1}$). These important references are almost not discussed, as Part 2.2 is kind of presented as an original contribution, whereas it is a small modification of very classical algorithms (whose original ideas trace back to Gower in 1966 [2]). It seems to me that this method is actually doing a simple MDS but in the embedding space defined by the kernel. To make these connections clearer, I think it is important to derive the loss that is minimized by the algorithm.\n\nAnother point is that there is no guarantee that the presented approach actually leads to a good estimator of the true kernel $K$ compared to classical maximum likehood estimators (as in Gaussian processes). It would be really interesting to prove that the estimator proposed is somehow consitent and is able to estimate the true latent positions. \n\nMoreover, a recurring claim in this article is the algorithm's speed, yet there is no analysis of algorithmic complexity. It doesn't seem that a solution based on eigenvalue decomposition is especially efficient, as it is of the order of $O(T^3)$ in computation. Concerning this point, for displaying runtimes, a logarithmic scale on the y-axis would provide a better idea of the differences between methods because, in Figure 2, they all appear roughly identical (except GPLVM). Additionally, for a fair comparison and to claim that \"IKD is able to provide competitive performance against some SOTA optimization-based methods but at a much faster running speed\" it would be necessary to calculate the speed of methods for varying $T$ and not just for a dimension $N$ that varies.\n\nFinally, the experiments don't seem to support the claim that \"As an eigen-decomposition-based method, IKD achieves more reasonable latent representations\" and \"IKD always obtains improved performance with an increasing observation dimensionality, and claims its absolute superiority when the observation dimensionality is very large.\" Indeed, the only metric proposed in the experiments is a classification metric with a $k$-NN on the representations. To evaluate the quality of an embedding and claim that \"IKD always obtains improved performance\" this doesn't seem sufficient to me. It's also important to compare methods with other metrics such as the silhouette score or trustworthiness. Additionally, Figure 6 is quite challenging to interpret, making it difficult to draw any conclusions. Variance is not presented, the curves are hard to read, and the significance of different $k$ values is unclear (the three lines are almost identical). Here, an aggregated table seems more appropriate, with a setting such as $k=1$ for example. Finally, the experimental setup for competing methods is not entirely clear, especially since t-SNE is sensitive to the perplexity parameter, and it's unclear how this parameter is chosen.\n\n[1] Euclidean Distance Matrices. Essential Theory, Algorithms and Applications. Ivan Dokmanic, Reza Parhizkar, Juri Ranieri and Martin Vetterli.\n\n[2] Gower, John C. Some distance properties of latent root and vector methods used in multivariate analysis, 1966.\n\n[3] Multidimensional Scaling, Sammon Mapping, and Isomap: Tutorial and Survey, Benyamin Ghojogh, Ali Ghodsi, Fakhri Karray, Mark Crowley.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This article presents a dimension reduction method based on estimating latent variables from a kernel derived from observations. The core idea is to see the empirical covaraince matrix as an estimation of a kernel matrix between latent variables (under Gaussian assumption) and to try to estimate this latent variables with eigenvalue decompostion. The original idea of considering the intermediate step of a Gaussian process for dimension reduction is intersting a seems to work reasonably well in practice.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Improve the experiments as detailed above + provide more clarity in establishing connections with previous standard spectral methods for dimension reduction.", "broader_impact_concerns": "No ethical implications of the work.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "5nVJ9PZQTQ", "reviewer": "Reviewer_9zQM", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths: The numerical results, especially for real data experiments are good, demonstrating the advantages of IKD, compared with PCA, KPCA, LE, Isomap, t-SNE, UMAP, GPLVM and VAE. \n\nWeaknesses:\n\n1. The paper is not written in a good shape. For example, there lacks several words to introduce the ``problem setting''. In particular, how to measure the quality of a dimensionality reduction is not presented. \n\n2. The authors argue that their approach is capable of tackling noisy data. However, in their simulations, the noise is too small, i.e. \\mathcal N(0,0.05^2).\n\n3. The motivation and algorithm are not written well. There are so many tedious but simple mathematical computations that are easy to derive.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper focuses on dimensionality reduction problems. The authors propose a novel nonlinear dimensionality reduction method, called as Inverse kernel decomposition (IKD), for this purpose. Several numerical simulations have been conducted to illustrute the effectiveness and outperformance of IKD.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "1. The paper should be rewritten, especially to highlight the motivation (why IKD performs well rather than exhibiting tedious mathematical expressions), algorithm itself (how the algorithm implementation, why the specific kernels are adopted) and the explanations of the algorithms (the effect of algorithmic parameters, the computation burden, etc..).\n\n2. Highlight the problem setting, especially demonstrating how to measure the quality of the mentioned approaches to show the advantage of the proposed IKD. \n\n3. Experiment: 1) The effect of noise, or the outperformance of IKD, should be added.\n                         2) The selection of algorithm parameters (or kernels) and the role of the algorithm parameters should be given.", "broader_impact_concerns": "NA", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "mBfLrLPqKB", "reviewer": "Reviewer_Mh8Z", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength:\n\n1. The proposed inverse kernel decomposition seems new to the best of my knowledge.\n\n2. The proposed method is simple and efficient. \n\n3. The empirical performance is promising.\n\nWeaknesses:\n\n1. Given so many eigen-decomposition-based dimensionality reduction method, the novelty of this paper seem rather limited. \n\n2. The motivation is not clearly stated. Why do the authors need to rely the dimensionality reduction method on the Gaussian process latent variable model? \n\n3. The writing is not always satisfactory. For example, there are typos like \"all pairs of rows in Z. I.e., ...\"", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper proposes inverse kernel decomposition (IKD), a dimensionality reduction method based on the Gaussian process latent variable model. The proposed method is simple and efficient. Promising experimental results are reported on synthetic and real data.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "It is suggested to revise the manuscript according to the above Weaknesses.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper introduces a novel nonlinear dimensionality reduction method named Inverse Kernel Decomposition (IKD), which relies on an eigen-decomposition of the sample covariance matrix of data and draws inspiration from Gaussian process latent variable models. Reviewers have acknowledged the novelty of the proposed method. **Two reviewers**  have strongly recommended clear acceptance following the rebuttal with the authors.

**One reviewer** commended the authors' efforts in enhancing the manuscript during the rebuttal; however, they expressed some conservatism regarding the current version. Specifically, they suggested two points for improvement:

1) Enhancing the readability of the experiments by ensuring clarity in the figures and explanations of the conclusions.
2) Providing more extensive discussions on the connections with prior research and existing spectral methods. While the authors mention Isomap, it appears that their approach may be viewed as a specific instance of kernel PCA with a particular choice of kernel.

Considering these suggestions, I recommend acceptance with minor revisions. I encourage the authors to address the aforementioned comments to further strengthen the paper.
Decision: Accept (poster)