Dynamics of the accelerated t-SNE
Authors: Kyoichi Iwasaki
Hideitsu Hino
Venue: Accepted by TMLR ()

Abstract: This paper investigates the dynamics of t-Stochastic Neighbor Embedding (t-SNE), a popular tool for visualizing complex datasets in exploratory data analysis, optimized by the Nesterov’s accelerated gradient method. Building on the foundational work that connects t-SNE with spectral clustering and dynamical systems, we extend the analysis to include accelerated dynamics which is not addressed in the previous work, revealing the emergence of Bessel and modified Bessel functions as a novel aspect of the algorithm’s behavior characterizing the temporal evolution of the accelerated t-SNE. Because the ordinary differential equation corresponding to the optimization process under consideration has a closed-form solution, by performing eigenvalue decomposition of the data’s adjacency matrix as a pre-processing step, we can obtain low-dimensional embeddings at any point in time without performing sequential optimization. This advancement not only enhances the practical utility of t-SNE but also contributes to a deeper understanding of its underlying dynamics.

Published in Transactions on Machine Learning Research (06/2025)
Dynamics of the accelerated t-SNE
Kyoichi Iwasaki
amaaskki.kyoicchi@ssii
The Groduate University for Advanced Studies, SOKENDAI
Hideitsu Hino
hino @ism. ac.jp
The Institute of Statisticnl Mathematics
Reviewed OpenReview: https: // openreview. net/ forum? 42=5
huttpe://ppenreviee/
Abstract
This paper investigates the dynamics of t-Stochastic Neighbor Embedding (t-SNE), a pop- ular tool for visualizing complex datasets in exploratory data analysis, optimized by the Nesterov's accelerated gradient method. Building on the foundational work that connects t-SNE with spectral clustering and dynamical systems, we extend the analysis to include accelerated dynamics which not addressed in the previous work, revealing the emergence of Bessel and modified Bessel functions as novel aspect of the algorithm's behavior charac- terizing the temporal evolution of the accelerated t-SNE. Because the ordinary differential equalion corresponding to the optimization process under consideration has closed-form solution, by performing eigenvalue decomposition of the data's adjacency matrix as a pre- processing step, we can obtain low-dimensional emboddings at any point in time without performing sequential optimization. This advancement not only enhances the practical util- ity of I-SNE but also contributes to deeper understanding of its underlying dynamics.
1 Introduction
Data visualization plays pivotal role in exploratory data analysis, enabling the comprehension of complex datascts by reducing them to cognitively managcable two threo-dimensional representations. Among various dimensionality reduction techniques, I-distributed stochastic neighbor embedding(L-SNE) (van der Maaten & Hinton 2008), descendant of the stochastic neighbor embedding (Hinton & Roweis 2002), has emerged as powerful tools for visualizing the overall structure data and understanding how clusters are formed across diverse fields. (-SNE continues to undergo technical and theoretical development, and its applications are also expanding. Tn (Zhu & Ting. 2022), the use of an isolation kernel instead of Gaussian kernel accelerates computation, while (Böhm al. 2023) integrates contrastive learning within computer vision rescarch, and (Skrodzki al. 2023) explores the concept of coarse embedding to expedite inference.
Furthermore, the concept has been extended into hyperbolic space using polar quadtree techniques,
demonstrated in related research by {Zhou & Sharpee 2021 Guo et al. 2022). Also, theoretical advancements have been made, such as the convergence study presented in (Jeong & Wu 2024). For general low-dinensional representation, the concept dynamical dimension reduction is introduced in Yoon & Osting 2022
In this paper, we focus on another aspoct of t-SNE, the acceleration of the optimization procodure from the viewpoint of dynamical system. In many software implementation of (-SNE (Szymański & Kajdanow- iez 2017 Krijthe 2015), the optimization procedure includes momentum term. We consider applying momentum-based methods such as the Momentum Method (MM) and the Nesterov's accelerated gradient (NAG) method (Nesterov) 1983) to t-SNE aiming to enhance the convergence behavior of the optimiza- tion dynamics in the early exaggeration (EE) stage. We analyze the impact of these modifications from a continuous- dynamical systems perspective.
We note that the t-SNE is composed of two stages: the EE stage and the embedding stage (van der Maaten & Hinton 2008). Some previous theoretical studies (Linderman & Steinerberger 2019 Cai & Ma 2022 have

---

Published in Transactions on Machine Learning Research (06/2025)
conducted analyses focusing on the EE stage. While these methods are often introduced to accelerate prac- tical optimization, our primary focus here is on understanding the dynamical behavior of the system rather than achieving computational speed-up. Since this stage itself possesses sufficiently rich mathematical structure, we will also focus on the EE stage, and the embedding stage is not modified in this study.
1.1 Contributions
We draw upon the foundational work of (Cai & Ma 2022): who established a connection between t-SNE clustering dynamics and spectral clustering under certain conditions for gradient descent (GD) method, framing the data clustering phenomenon within a dynamical systems perspective. While the use of MM or NAG can potentially accelerate the optimization in the EE stage, it is not obvious whether the resulting embeddings still preserve cluster structures in the same way as standard GD. Our contribution extends this analysis by demonstrating that the incorporation of MM and NAG can accelerate the EE dynamics; in particular, show that NAG introduces mathematical entities such the Bessel and modified Bessel functions into the system. Through rigorous theoretical analysis and empirical validation, we show that our approach offers a further insights into the dynamical properties of t-SNE.
Our main contributions are summarized as follows:
We explore an approach of dynamical system analysis of t-SNE including a momentum term, and further extended to the NAG. We find that in the former case (MM), even with the same linear ODE, convergence is theoretically faster [Proposition 5.4 and in the case of the NAG, it leads to the derivation of related ODEs as the first kind Bessel functions or modified Bessel functions. [Theorem 6.4 As these claims suggest, our approach explicitly involves eigenvalues and eigenvectors. When the data size is n, solving the eigenvalue problem requires computational complexity O(n"), which is greater than the original On²²) complexity. Therefore, in this paper, we limit our numerical experiments to datasets of relatively small volume.
We derive these ODEs and their relations, organizing them into several theorems, propositions and lemmas, where all proofs are provided the Appendix. As matter of fact, deriving parallel results already given without momentum in (Cai & Ma 2022) to those with momentum term are technically not straight forward.
Since the closed-form solution of the ODE is available, embeddings at any point in time can be easily obtained. [Proposition 5.4 Theorem 6.4
We find that insights related to EE stage in (Cai & Ma 2022 are replicated in our method with acceleration, leading to the conclusion that the optimization algorithm converges in less time than when using GD. Since the EE stage does not have an explicitly defined objective function with the acceleration parameter a > 1, it is challenging to monitor the function value and determine an ap- propriate stopping point during optimization. As a secondary benefit of our continuous formulation via ODEs and their closed-form solutions, we can define a principled stopping criterion based on a threshold applied to the Accumulated Residual Ratio (ARR) of major and minor eigenvalue contri- butions. (27). (2 This approach enables the identification of a suitable termination time within the continuous franework. Such mechanism is particularly advantagcous when t-SNE is employed as part of a larger system or pipeline, where an automatic and interpretable stopping condition is beneficial.
We theoretically demonstrate that, under certain conditions, iterative computations for both MM and NAG can be replaced by their continuous counterparts by adopting this continuous relaxation approach using ODEs. [Proposition 5.2 6.2 It provides thcoretical justification for results that had previously been obtained heuristically through ilerative methods. We can confidently rely on the original discrete optimization algorithms (MM and NAG), as they retain the essential dynamics while being more computationally efficient. However, since solving the ODEs typically requires cigenvalue decompositions of large matrices, the practical advantage of using the continuous formulation is limited.
2

---

Published in Transactions on Machine Learning Research (06/2025)
By interpreting the spectral decompasition of the dynamics in terms of Lyapunov exponents, offer novel theoretical perspective on the formation of clusters in t-SNE. [Proposition H.1 H.2 Components associated with positive Lyapunov exponents dominate the long-term behavior and effectively determine the clustering structure, with the corresponding cigenvoctors encoding cluster assignments. ((2)) In contrast, components with negative Lyapunov exponents decay over time and gradually vanish from the embedding. Notably, the case of NAG, these stable components do not simply decay exponentially but exhibit non-monotonic, oscillatory behavior governed by Bessel functions.
1.2 Key take aways for t-SNE users
We investigate both the theoretical and practical aspects of MM and NAG, based on the framework of GD. Our study builds upon the theoretical foundation established by (Cai & Ma 2022 which demonstrated that the EE stage of the t-SNE algorithm can be approximated by a first-order ODE. The solution of this ODE can be expressed explicitly terms of the adjacency matrix's cigenvalues and cigenvectors, and providing a theoretical justification for why I-SNE reveals the data's cluster struclure.
In our paper, extend this line of reasoning by investigating whether similar theoretical framework be applied to other optimization schemes such as MM and NAG. Specifically, we ask: Can the dynamics of MM and NAG during the EE stage also be described using ODEs, and if so, do the same assurances hold?
The answer to this question "yes": we theoretically demonstrate that optimizing t-SNE with either MM or NAG likewise yields the emergence of clear cluster structures. This result implics that, although the computational cost of optimization was previously concern for t-SNE, one now employ accelerated gradient methods with confidence. Furthermore, our work proposes novel stopping criterion for the ODE- based formulation, point which is not explicitly addressed in (Cai & Ma 2022)'s 's original work. This contributes both to the theoretical understanding and practical usability of accelerated variants of t-SNE.
Next, we present the practitioner-oriented takeaway. We derive the corresponding ODEs and their closed- form solutions by continuously relaxing the step-by-step ilerative process. This relaxation enables, particu- larly in the EE stage, the approximation of intermediate states arbitrary time points by evaluating the ODE solution. In addition, we propose a stopping criterion, referred to as the ARR metric defined in (28 which provides a principled method for determining an appropriate termination time. To apply our approach in practice, two major assumptions must be taken into account:
When the data size n is sufficiently small, the cigenvalue problem required by the method method-typically of computational complexity O(x") can be solved within a practical amount of time.
The assumptions stated in the theoretical takeaway-nande (II), (TLM),, (T2.M) for the MM method, and (TI), (TINN), (T2.N) for NAG are satisfied.
Provided these conditions hold, the following approach becomes applicable:
Preparation: Compute the adjacency matrix and its eigenvalues.
EE stage: Compute the ARR-based criterion and determine the stopping time.
Embedding stage: Obtain the low-dimensional embedding by evaluating the formula (3) with mlk+1) 0, where acceleration unnecessary in order to obtain stable results in the embedding
stage.
2
Related work
Dealing with t-SNE, a method of dimensionality reduction, as a dynamical system requires a background in both the methodology of dimensionality reduction and analysis of algorithms as dynamical systems. Moreover, the unique properties of t-SNE and the development of related research are also pertinent. We concisely summarize the related studies for each aspect..
3

---

Published in Transactions on Machine Learning Research (06/2025)
2.1 Dimensional reduction methodology
In dimensionality reduction, numerous methods have been established to transform high-dimensional data into a more managcable, lower-dimensional representation, aiding in clucidation of the data structure. Prin- cipal component analysis (PCA) (Joliile 1986) serves as a foundational linear method, identifying principal axes of maximum variance. Isometric mapping (Isomap) (Tenenbaum et al. 2000) focuses preserving the geometric distances on the manifold, effectively maintaining the spatial relationships among data points. Locally linear embedding (LLE) (Roweis & Saul 2000) isolates local linear relationships by reconstructing points from their immediate neighbors, reflecting the intrinsic geometry of the manifold. The Laplacian eigenmap technique (Belkin Niyogi, 2001) constructs a graph Laplacian and leverages its eigenfunctions to achieve a lower-dimensionall projection that reveals the manifold's local characteristics.
By further developing graph-based methods such Isomap and LLF., SNE (Hinton & Roweis 2002) and t-SNE (van der Maaten & Hinton 2008) have been proposed as methods that learn embeddings directly and have become very popular. UMAP (Melunes & Healy 2018) and PHATE (Moon et al. 2019) have been developed as powerful alternatives to the t-SNE algorithm, and recently, the relationships between t-SNE
and UMAP have also been discussed Damrich et al. 2023
2.2 Gradient methods to differential equation
The approach of deriving and analyzing ODEs as continuous limits of ilerative optimization algorithms is widely adopted. This includes attempts to extend acceleration techniques beyond NAG (Wibisono et al. 2016), employing Runge-Kutta integrators (Zhang et al., 2018) and symplectic integrators (Goto & Hino, 2025), and formulating continuous-time models using Lyapunov analysis and tools from stochastic calcu- lus (Orvieto & Lucchi 2019). These studies exemplify the breadth of research exploring different aspects of continuous approaches to optimization. There are more geometrically sophisticated methods, such as (De- fazio 2019). which addresses NAG on Riemannian manifold, and (Wilson et al. 2021), which performs analysis using Lyapunov functions. In deriving the ODEs related to the NAG method, we refer to several results from [Su et al. 2016), where the emergence of Bessel functions for a related ODE is established.
2.3 SNE algorithm and recent developments
In addition to the studics introduced in the introduction, variants of t-SNE and relationships with other methods are also being actively researched (Kobak et al. 2019 Chen et al. 2015 Pezzolti et al.] 2017 Fujiwara et al., 2020). A generic formulation of embedding algorithms that includes SNE and other existing algorithms are studied, clucidating their relation with spectral mothods and graph Laplacians (Carreira-
Perpiñán 2010
Vladymyrovv & Carreira-Perpiñann 2012).
There are several theoretical attempts, such as [Linderman & Steinerberger 2019). which compares spectral clustering and t-SNE and (Arora et al. 2018), which states that t-SNE works well under specific conditions called 7-spherical and y-separated in EE stage. The paper (Linderman & Steinerberger, 2022 presents theoretical and experimental results, focusing on sufficient conditions for the parameters for acceleration, and h for step size under the assumption of GD. It also discusses considerations regarding disjoint clusters and the independence of the results from initial values. Note that such results include both EE and subsequent embedding stages.
Additionally, recent studies have addressed aspects related to implicit regularization and scaling effects in clustering algorithms. For instance, (Aufinger & Fletcher 2023) discusses implicit regularization mecha- nisms, as highlighted in Proposition 6.5 which connects closely to early stopping in clustering optimization procosses. Additionally, (Murray & Pickarski 2024) explores the scaling effects arising from the gradual decrease of Pij, which is the element of adjacency matrix commonly used in t-SNE is the symmetric proba- bility matrix P = (p.y), as n + cc. This phenomenon appears to relate to Lemma 5.5 or 6.3 particularly the condition concerning the sign of (LR)L(P)) - ~ where AR()(P) stands for Etth eigenvallu of the unnormalized Laplacian matrix of P. These findings provide critical insights into how regularization and scaling behaviors influence clustering performance and stability under large-scale settings.
4

---

Published in Transactions on Machine Learning Research (06/2025)
2.4 Acceleration methods
It is widely known that drawback of L-SNE is its high computational cost, and many attempls have been made to improve its computational efficiency. The original paper (van der Maaten & Hinton 2008) uses both early exaggeration effect by adopting acceleration parameter a >> 1 and momentum methods. The Barnes- Hut tree is a well-known approach for accelerating the t-SNE optimization procedure (van der Maaten 2013 2014). The n-body force calculations is utilized in (Yang et al.; 2013 Vladymyrov & Carreira-Perppiann 2014]. (Skrodzki al. (2024) applies Barnes-Hut approximation idea into hyperbolic disc to accelerate low- dimensional embedding. We note that is experimentally shown in (Lambert, et. al. 2023 that NAG can be an alternative to EE
Although the term "acceleration" often suggests computational speed-up, we clarify that our current frame- work focuses on theoretical acceleration in the continuous-time dynamics, rather than practical runtime efficiency, In fact, due to the reliance on full eigendecomposition of the graph Laplacian, the overall com- putational complexity remains O(n²) in our current implementation. Nonctheless, we believe that this spectral formulation offers valuable analytical insight, and that future implementations using approximate eigensolvers (e.g., Lanczos algorithm) may help bridge the gap between theoretical structure and practical scalability.
For example, the oscillatory behavior revealed through Bessel-type dynamics is directly captured in the ARR metric we propose, which provides an interpretable and data-driven criterion for early stopping. This ana- lytical formulation could inspire new techniques for adaptive exaggeration scheduling or spectrally informed initialization, potentially improving both the interpretability and efficiency of t-SNE variants in largo-scale applications.
3 Original t-SNE algorithm
We explain the original t-SNE algorithm fvan der Maaten & Hinton 2008
For a vector a = (at)isisn E R", define its I' norm by ||a||p = (E la P).P. For a matrix A = (a:y) E Raxn and its spectral norm is
||A|| x ||Ax||2 We also introduce sets of orthogonal matrices O(n, k) = =V R=*k, V TTy V [
and (nn) = (m,, n). We use the notation [n] [n]=60 (1,,,, ..2....,,,, n} for Z The diameter of set S(C R*) is
diam(S) (S) s. || yl|2- For symmetric matrix A = (uyy)ssssnn E Bu we define degree operator D R - R D(A)= Hiag(E: ===1 da.,.... Ex= aia), and Laplacian operator L Rexn - Baxn by L(A)- D(A) A. Also, we set 1, = 1) E R" and H (1, - I.) E g"x*.
Let (X.) C Rd be the dimdinnnsional data points. For the pairs of data points {(x,, Xy)hsttjnn we define a symmetric matrix P = (Pij) E Raxn with pa 0 for Vi [n] and pij (Pats + (2n) for # j, where
exp(- X - X, 1
Pfl Eie( 1|
X 1B/2T2)
(1)
Here, T; > 0 is the bandwidth parameter. The aim of t-SNE is to achieve a low-dimensional representation fyr (ee)e S Rd, di < d, of (Xiheeo)) Typically, dj is two or three because t-SNE is predominantly used for visualization, and in this paper, consider d without loss of generality. Consider symmetric matrix Q (gr))ssiggn where qi = 0 for Vi [n] and
(1 + 1b:
4i Σ (1.2,....)) (1 + 11 y. 111
(2)
for # j. We get the low dimensional representation (MDee by minimizing the Kullback- Leibler (KL.) divergence Kullback & Leibler 1951) between matrices P and as (v.,..., Uy) arg min Σi,jeln) (æj Pij log PS, 933
As described in (Cai & Ma 2022), the t-SNE algorithm typically implemented in two distinct phases: an EE stage followed by an embedding stage. While previous studies have explored accelerating the EE stage by introducing parameter
into the GD framework, this paper extends the concept further
5

---

Published in Transactions on Machine Learning Research (06/2025)
by investigating whether analogous acceleration strategies can be systematically applied to MM and NAG. Although the explicit form of the objective function for the case with > 1 is not known, proceed by assumiin the following formulation in line with the original t-SNE work (van der Maaten & Hinton 2008 and Caa Ma 2022). The following update rule aligns to (Cai & Ma. 2022), which states
(1+1)
(k) - y/ +
E
(b) )5 (ar + (k++)) (uf yi ) yi
Vi
ISWJEn 188
where h R+ represents the step size parameter, originally denoted as 4h in ovan der Maaten & Hinton 2008), and m(k-1) R is momentum parameter. The value 0, known the exaggeration parameter, plays a crucial role in the gradient coefficient $/ (a) derived in Appendix A and is given by
(k)
(a) = 1 11 Dj 12
aPi 4i
E R
(4)
(9) where (0, - ( (a) E R², and we define adjacency matrix s (5 (a)) The
algorithm starts with initialization vi (0) y -11 vi [ni]. Although the momentum term m 1 (k)
yi (R-1) to accolerate the updating algorithm was not discussed for simplicity in (Cai Ma 2022), our aim to discuss how we may accelerate the convergence of the algorithm and derive a dynamical system corresponding to accelerated updating methods.
4 Asymptotic behavior of the update equation
The points {u, (8) hejmi C IR2 are represented in two dimensions, but they are reorganized into n vectors according to the first and second coordinates as y R., E [2]. This identification and the following Theorem 4.1 enable to make sure that the adjacency matrices s(k) is considered fixed matrix HH-H,, For 1 E [2) Eq. is reformulated
(k+1) = [I, - hL((S()))) (M) + m (x - --(111) -xx-)) yi
(5)
where In E R is the identity matrix and y/ah) E R" is the l-th coordinates of {v. (2) he(m) The following theorem shows that the matrix !! = ( (a))a.se(x) can be approximately simplified. This approxima- tion, combined with a continuous relaxation with Hn, enables the derivation of ODEs from the recurrence relation (3), making it one of the core statements of this paper.
Theorem 4.1. (Asymptotic Graphical interpretation, Theorem in (Cai & Mo 2022 Let P (pt))
be symmetric matriz defined in Eq. [1] and denote (T()) == (diam ((v, with + j, and each k 2 1 such that s7t < 1, we have
(2)
))s.sn))? Then, for any i, E [n]
(s (a) apu + n(n 1)
1
S apij) + n(n 11(1111111 17 ())
217(
(6)
Then for each k 2 1, as long as (8) and satisfy 1k)
PP||
> "|l we have
||S (aP H
lim
=0.
(7)
||P
H.ll
With this statement, Eq. (5 admits an approximation
[n hL(aP H,,)y + m ( y/ ((k-1))).
(8)
We consider the global behavior on {y]" (k)) l [2] under the following initialization and condition of param- eters H
The condition (TI) comes from 'Initiaiiaation, and (TI) describes 'Trajectory' as + 2o.

---

Published in Transactions on Machine Learning Research (06/2025)
(I1) (v) maxiepa lyy (0) Ioo. [ly]" O(1) 00.
Wy/hemn hem {ur (1) Jeen satisfy mingetal {1ly) (0) 2|2 lly? Ay/"11a [12] 0 and
(T)) The parameters (a, h k) satisfy k(nho |acbao||||ll P|l| + h//) 0(1) as n 000
Condition (II) states that the initialization { should not be simply all zeros or unbounded, and (T1)
says that the cumulative differences of hL( and hL(aP H.) are limited.
The update formula (8) is three-term recurrence relation and two initial vectors y and yt need to be set. With the momentum term, the vectors y: y (++1) stays within limited area:
Proposition 4.2. (Localization) Under the condition (I1) and (T1), we have
diam(Ly, ( S (ax |I b0 11y/")1>>))
(9)
[e[2]
for some universal conslant C 0.
The proposition 4.2 ensures that the data points (y are globally bounded and the condition (II) leads to (*)) < and both the theorem 4.1 and the approximate Eq. (8) holds. This result had been derived in (Cai
& Ma 2022) without the momentum term; however, deriving this result when considering the momentum term is not trivial and requires an evaluation of the difference in solutions at two consecutive points t - 1 and originating from the momentum lerm.
5 Gradient flow with constant momentum coefficient
We assume that m (b+1) = m E (0, 1), which is commonly assumed in existing works KKvachki & Stuart,
2021 He et al. 2023 Hao et al. 2021 Sutskever et al. 2013 Rashidi et al. 2020 Let (5) Jx20 be computed by the update equation
(K)
y/ (k+1)
[I, - ML(oP H²) + m(s y(k-11)
(10)
The following lemma derives the ODE that holds for MM the tontimuous-time limit as - Lemma 5.1. With the assumption (kk) Yi(kh), kh and the continuous timit h + 0, the update rule 110 is reduced to the ODE with initial value y(0) = Y.(
Y(() 1-m L(oP H,().
(11)
Furthermore, the following proposition provides a theoretical guarantee that 8/ and (ich)) become sulli- ciently close as h + 0.
Proposition 5.2. (Gradient flow with constant momentum) Under the condition (I1), for 1 E [2], Zzo, and some positive constants Ci, C2, we have
(M) Y2(kh)22 (2(0))22
2khmC1
k h _) 2 (1- m) ([(PP-HH)1P..
[[(P H.)ll
(12)
m
This proposition justifies that the approximation y Y, works under the conditions (I1), (T1). To
better understand the solution of the derived ODE that the matrix P is symmetric, and its Laplacian P is likewise. With the eigendccomposition of L(P) expressed L(P) VTAV, we have
we consider eigendccomposition. Remember

---

Published in Transactions on Machine Learning Research (06/2025)
Lemma 5.3. The matric L(aP- Hn) is symmetric, and we have eigendecomposition:
L(aP -
(13)
where Σ = diag(o:))e(n)) A = (A))g)))) 0= \1 A2 ... S Xn are eigenvalue and V E O(n), which are the same eigenvectors of L(P). Here, the relation between and is
σ1 = A 0, 0,0, a
72 (2 S n).
(14)
With the Lemma 5.1 we achieve the explicit expression
Yi( = exp 1 m
L(aP LLa
-
yf
(15)
Combining the eigendecompasition and Eq.. 15 yields the following proposition:
Proposition 5.4. (Solution path constant momentum) For I [2], the first order linear ODE
11 with
initial value (0))) = Yi and momentum term m E (0,1) has the unique solution:
(1) = (u yi mu Σ exp (-1-
1 I m ad: n ) (n y ) ).
1
(16)
Propositions 5.2 and
5.4 exlend (21) or Proposition 8 in Cai & Ma 2022 with constant momentum
term, as substituting m in Eq. 16 recovers the original result. The effect of momentum coefficients is discussed in Appendix L.7
5.1 Well-conditioned matrix
The behavior of exp (-r-m (a): - - ) for t >> depends on the sign of aλ; - If negative, the corresponding components are amplified: positive, they are suppressed. It implics that
lim Yi(t) E span((u 100
(17)
If the data {Xihsenn are well clustered and bandwidths 7, properly chosen, (Balakrishnan al. 2011 demonstrates that P is well approximated. We assume the weighted graph of well-conditioned matrix P as detailed in Appendix E has R 2 2 connected components.
5.2 Spectral convergence with constant momentum coefficient
Condition((1) ensures Theorem 4.1 holds. We assume stronger, analogous condition (TI.M), along with (T2.M), which governs cigenvalues and clustering behavior in section 5.1
(TL.M) The parameters (a, h, satisfy a >> [nAn+1 (L(P)) and + = o(n) as n oc.
I-m
(T2.M) There exists symmetric and well-conditioned matrix P E R such that (L=1)( I-m [(PP - = (11) as n oc.
max max(( ((1-=))) [L(P" - P)I) and to
In the context of visualization and clustering, the following lemina Jundamental identifying which components remain dominant in the eigendecomposition representation. Lemma 5.5. Under the conditions (T1.M), (T2.M) and > 1, we have
1 aAR(L(P)) n 00, 50, adR+1 (L()) n
1
0.
(18)
Especially, the eigencalues {a.) C I consist of R-negative numbers or zeros (0 = 01, 02,..,,,00 posétive numbers (ont1....,ow)
RR), and n R
We denote i-th eigenvalue as X(A); we emphasize original matrix
8

---

Published in Transactions on Machine Learning Research (06/2025)
Note that this proposition is not explicitly stated in (Cai & Ma 2022). but it implics that the conditions (TI.M) and (T2.M). This, in turn, means that the same R eigenvalues and eigenvectors in the GD case with m = 0 remain in the visualization, thereby providing theoretical guarantee that MM serves as an accoleration of GD.
18 holds under
6 Gradient flow with NAG
Another accelerated gradient method is NAG NNeserrov 1983 with ) and w/ y1)) under our situation:
which takes the following form: starting
(N) hL(PP Hn)w/ (8) 81 wi
(k+1)
(19)
w (6) (k)
yi
k - (k)
h: + ( y((-111
(20)
Similar to Lemnaa 5.1 have another linear ODE with NAG.
Lemma 6.1. With the assumption 81ah Yi(k/)) for smooth curve Y (t for and 1
[2], and
putting = /h, we have the ODE corresponding to the update rule of NAG in Eqs. 19 20 as
Y((t) + Y + Y,()) 0.
(21)
The following proposition provides a theoretical guarantee that y (k) Proposition 5.2
and Yi(k/)) holds for NAG, just
Proposition 6.2. As the step size h 0, the update rule that for all fized T > 0,
19
20 is reduced to the ODE 21 in the sense
lim max
ly (k)
Yi(kvh))l h)1l2 =
(22)
h h-0 0ocks
We aim to demonstrate that this situation would be reasonable under the following conditions:
(TI.N) The parameters (a, h) satisfy a > [nAR=1(L(P)] and = o(n2) as n - ac.
(T2.N) There is a symmetric and well-conditioned matrix P E R² such that AR+1(L(P!) > max((faa)=1 [L(P*
|LL(P PI))] and LP - P)]l = o(1) as n + oc.
Similar to Lemma 5.5 we have the following lemma:
Lemma 6.3. Under the conditions (TI.N), (T2.N) and n >> 1, we have
aAn(L(P)) - n 50, aAR+1(((P)) n - > 0.
(23)
Especially, the entire eigenvalues (o.....,o,) consist of R-negative numbers or zeros (0 = 0i, O2...... ,n)
and n - positive numbers (IR+I)...,,w,, on))
In the case of NAG, this proposition guarantees that, under the conditions (TLN) and (T2.N), the same R eigenvalues and eigenvectors emphasized in visualization and clustering, and that the number R matches that in both GD and MM as Lemma 5.5
As in Lemma 5.3 we have another dlosed-form expression by solving linear ODE
Theorem 6.4. Under the conditions (T1NN) (T2NN), partition the eigenvectors as V (U,U1), where U E On, R), U. E O(n R) Then, we have
Yi(1) = UT (t) + U (
(24)
T is orthogonal complement, of

---

Published in Transactions on Machine Learning Research (06/2025)
where
2 Ti(1) = diag 1, W-o 11(1-110 01)
1111 oR) .I2(1) diag 2 (1((/ott)))...... (WORI
2 2n(((von) Ivon
W--OR
and the symbol I1(-) denotes the modified Bessel function of the first kind, while J1(+) represents the Bessel (0) function of the first kind. Also, Yi,I. (yi, . 1 .i.R RR, (0) the inãtial values.
NiR++im = (yiR " ((n) g=-R are
The explicit expression 24 is the solution path under the NAG for the t-SNE algorithm. In the Appendix L1 we consider the qualitative differences in the solutions of ODEs corresponding to GD, MM, and NAG by plotting the solutions.
Lemmas 5.5 and 6.3 demonstrate that eigenvalues and eigenvectors for which ad; - - S 0 are emphasized,
highlighting key features for spectral clustering. This phenomenon, where large eigenvalues X decrease with increasing is referred to implicit regularization. The subsequent proposition asserts that Yi(t) converges towards the subspaces spanned by the first R cigenvectors, which form the Laplacian mull space commonly utilized in spectral clustering, thus embodying implicit regularization with the inclusion of momentum terms. Furthermore, the assumptions regarding t such as t o(n), t o(n ta ( - = o(1) or ea Pa||L(P P)ll = o(1) imply necessity for early stopping to avoid overshooting.
Proposition 6.5. (Implicit regularization, clustering and early stopping) Under the conditions (I1), (TI.M)
and (T2.M), let (t))) be the function defined
16. Uo On, R) be such that the columns span the null
space of P Then we have
Yi(1) U. Y,(Oll = 0, LE [2].
lim 410000
(25)
Y
Also, for a permutation matriz E R"x,, we have
Y,(() -2ull2 0, LE [2]. (2(00)12
lim R-000
(26)
(0) / /n, for [R], the wumber where z (211.-.,, 211, 2125 212:---. 2UR....,2.R)) = E R and zir = 0 Vi
of By is nyy i.e. the ramberr of nodes in the r-th connected component. Similarty, under conditions (11), (TI.N) and (T2.N), let (i) be the function defined in (24) Uo E On, Rt be that the columns span the null
space of p Then we have the same formulas
25 266.
This proposition shows that in both cases, Yi(I) serves as eigenvectors defining the Laplacian mull space, producing results equivalent to spectral clustering. Despite utilizing the same data points and the same adjacency matrix P with corresponding eigenvalues and eigenvectors, the increase in t presupposes distinct formulations: the MM adheres to & = o(n) and to Lo - P)ll = o(1) as n - co, whereas the NAG relies o22 and Pa]]L(( P o(1) as + co. is consistent with the general observation that the NAG converges more rapidly than the MM. Also, Lyapunov exponents are discussed in Appendix H
7
Stop timing
While commonly used metrics such Trustworthiness (Venna & Kaski 2001 and Continuity (Kaski et al. 2003) are effective in evaluating the quality of low-dimensional embeddings after dimensionality reduction, they are not suitable for determining the evaluation timing in the context of ODEs. To address this limitation, we propose novel metric, Average Residual Ratio(ARR). As will be confirmed in Section 8.2 and [.6 the comparison between the proposed ARR and the Trustworthiness values of the embedding indicates well- aligned correspondence. The term 'stop time" is used informally to refer to a general criterion for deciding when to stop the algorithm, and should not. be confused with the technical notion of stopping time in stochastic analysis. With Eqs. (16). (24)), we express Y,(t) = [(t)=>ın, (t) with the time-variant coefficient (t)(t) R. To evaluate the positive contribution of ct,((t), define (t)t) = Eiery (c),(t)) With Eq. (14) and
in the case of MM, (1) = EieØ exp (-11g)
taz I-m
|u yi where we set m = for GD. By contrast, in the
10

---

Published in Transactions on Machine Learning Research (06/2025)
case of NAG,
2yL.
(0)
1 11 [ -oi).
(i= L...., R)
W--0
a1(1) =
1E|21
(27)
2y (0)
J1 (1
(i= R + 1,...,
Woo
GD (k-5,h-10,t-501
GD (k-22.2-10..-22)
GD (0-30.h-100--300
GD (K-70,h-10:t-700
0.08
0.06
0.02
0.00
-0.06
-0.08
-0.05
0.00
0.05
0.04
-0.02
0.00
0.02
0.04
-0.02
0.00
0.02
0000
-0.01
0.00
0.01
0.02
MM (k=5/h=10,==501
MM 1k=
MM (k=30.h=101=3300
MM (k=70,h=10,:=700)
-0.04
-0006
0.08
-0.05
0.00
0.05
-0.02
0.00
0.02
-0002
-0.01 0.00
0.01
0.02
-0.02-0.01 0.00 0.01 0.02
NAG (k-5,h-10,t-155
NAG (k- (k-22.h-10,t-69..
NAG (k-30,h-10.t-94.87
NAG (k=70)h-10,tt222 221.36)
Values smurf
0.10
neptune normal
0.05
satan
0.00
-0.10 0.005 0.00 0.05 0.10
0.00
-0.0
0.00
-0.5
Figure 1: t-SNE visualization of n = 500 KDD Cup 1999 dataset using three ODEs: GD, MM and NAG.
To determine the stopping time, we propose the Average Residual Ratio(ARR), which divides the equation
into main terms (i 1,.., R) and residuals (i R + 1,.., n), based on the average contribution to clustering:
セ- EiRn az(l)
ARR(() =
(28
R EE ai(t) + n- Ei= (1)))
B
11

---

Published in Transactions on Machine Learning Research (06/2025)
Based on the Lemmas 5.5
6.3
and the observation of Figures 4a 4b in Appendix, it is expected that ax(t)
increase monotonically for i 1,... R, and decreasing or gradually decreasing for i R + 1 1...,,
We propose a stopping criterion based on a threshold (e.g., 0.01) and demonstrate its effectiveness in subsec- tion 8.2 While t-SNE unsupervised, this heuristic is similar to those used for convergence in algorithms like GD MM/NAG. Our analysis shows that components with negative eigenvalues form clusters, while those with positive eigenvalues don't contribute the final embedding. Thus, using the ratio of these components convergence criterion for the ODE is reasonable.
7.1 Lyapunov exponents
The ARR used to determine the stop timing was based on the time parameter 4, but considering the solution of the ODE dynamical system reveals an intriguing relationship.
The Lyapunov exponents for MM and NAG are given by H.122 R eigenvalues are non-positive. Let the Lyapunov exponent for MM be ALgep, MM, and that for NAG be
H. R.123
Lemmas 5.5
6.3
show that the first
AAG.,,NAG, When momentum parameter m does not get close to 1, it is observed that
ALyep.MM. = 1- m (ad: n 1 <
<v aad= <v-a
n-1 ALgap,, for i [R].
(29)
This is guaranteed by the fact that the function of satisfies 0 x < 1 and x < VT, and larger Lyapunov exponent implics faster cluster convergence. This is consistent with the commonly stated notion that NAG achieves faster convergence optimization method.
8 Numerical experiments
We mumerically evaluate our results by comparing the t-SNE algorithms with and without acceleration and their ODE counterparts with real-world datasets: KDDeup1999 and MNIST Results on other data are shown in Appendix [ In this paper, we set perplexity = 30, as in Cai Ma (0022) In all experiments in this section, all initial embedding vectors were generated randomly. The discussion on the methods of initialization in t-SNE and the impact of different initializations, along with experimental results, are presented and discussed in Appendices and 1.5
8.1 Experiment 1: Dynamic behavior of clustering over time
The KDDCup1999 dataset (Stolfo & Chan 1999 Stolfo et al., 1998 Tavallace et al. 2009 includes wide variety of int simulated in a military network environment, and comprises 42 columns with their labels. We extracted 100 samples of data for each of the 5 labels (fsmurf', 'neptume', 'normal', bback' 'satan') from the dataset, totaling 500 samples.
For the solution of ODEs corresponding to GD, MM and NAG, we substitute concrete value of the time to obtain two dimensional map of the original points in 42 dimensional space in Fig. [ where the top, middle and bottom rows correspond same iteration count k, and we move to the right, k increases, with the corresponding time t increasing accordingly. The first row corresponds to GD, the second to MM, and the third to NAG. Comparing GD and MM, as described in Eq. 16). the effect of the momentum term m causes the results to appear though time /. has been fast-forwarded, even for the same time point. As for NAG, it can be observed that cluster formation occurs earlier in time t. Specifically, while the second column MM corresponds to t = 220, the fourth column of NAG at = 221.36 already shows that cluster formation is nearly complete. We state the assumption in Lemma 5. about GD and MM, that is we identify the variables as l kh, where the variables k is iteration number of the original iterative algorithms, h step size and is time variable. For NAG, the identification can be t = k as in Lemma 6.1
results of GD, MM and NAG, respectively. The column is aligned with the
4The experiments 16GB memory, using Python.
conducted
laptop PC with 12th Gen Intel(R) Corc(TM) 17-1255U procossor, 500GB storage,
á Sperplexity is typically between and 50, as mentioned van der Maaten Hinton
2008
12

---

Published in Transactions on Machine Learning Research (06/2025)
Let's pick up the first MM case (t 50) and the second NAG case (t 69.37). The variables are formally different, but for MM, due to the presence of the momentum term m = 0.5, it effectively becomes t = 69 50 - (See the formula in Proposition 5.4 where we adopt that MM can be understood as fast-forwarding GD with respoct to the variable t.). Although they have practically similar variable t, the results are drastically different. The three plots in the rightmost column are all elongated.
Furthermore, in NAG, the final plot on the rightmost panel corresponds to = 221.36, while in the MM, the second plot from the left corresponds to t = 314.29. The application conditions for the NAG are either t = 022 in (TINN and t² |P P PI|| o(1) condition in (T2.N). By contrast, the application conditions for the MM are either I-m = o(n) in (TI.M) and the condition I-m ||P - PIl o(1 in (T2.M). This generally suggests that the NAG can only be applied for relatively smaller values compared to the MM.
8.2 Experiment 2: Stopping criterion
Although we have conducted experiments to demonstrate that the accelcrated method can indoed achieve embeddings in fewer time steps in the previous section, in our analysis, the embeddings are explicitly obtained as solutions to differential equations. Therefore, while this research originated from accelerating t-SNE, it no longer makes sense to compare it with other methods from the perspective of computational efficiency, We demonstrate that by setting threshold to ARR (28), we obtain reasonable embedding without using an ilerative optimization algorithm with popular MNIST dataset, which contains grayscale images of handwritten digits, and we focus on 400 images with 4 labels("2", '4', '6', '8') totaling 1600 samples, where each image contains 28 x 28 = 784 pixels. At the same time, we provide comparison with standard Trustworthiness metric (Venna & Kaski 2001), which is commonly used in dimensionality reduction methods such as t-SNE.
Figure 2 shows the transition of ARR 28 and Trustworthiness with 10 nearest neighbors over the time (left) and the embedding results for each method: GD, MM, and NAG after early exaggeration stage (right-top) and embedding stage (right-bottom).
As increases, ARR decreases, with NAG being the fastest, followed by MM, and GD being the slowest. The ARR crosses the 0.01 threshold at = 2700 for GD, t = 1400 for MM, and t = 280.0 for NAG. The right panels show that the embeddings at these times provide reasonable clustering. As noted section 7.1. cluster formation occurs fastest with NAG, followed by MM and GD. This can be interpreted by analyzing the Lyapunov exponents of the underlying ODE. When comparing ARR and Trustworthiness, ARR decreases, Trustworthiness tends to increase, showing general inverse correlation. The correlation cocfficient was approximately -0.77. More details are reported in Appendix L.6.
As shown in the original paper (van der Maaten & Hinton, 2008) and the subsequent study in (Cai & Ma 2022), the t-SNE algorithm typically visualizes clustering results by performing an embedding stage Tollowing the early exaggeration stage, which intend to facilitate cluster formation. The visualization in the lower right of Figure 2 follows this convention. It important to note that the embedding stage is solely for visualization purposes. Therefore, no acceleration methods such as MM or NAG were applied. Instead, GD was iteratively used after the early exaggeration stage.
9 Conclusion
In this paper, we derive linear ODE and Bessel differential equation corresponding to the MM and NAG in optimizing t-SNE with explicit solutions. In addition, akin to results obtained through GD, perform- ing an cigenvalue decomposition confirms that implicit regularization, specifically clustering around small eigenvalues near zero, which is crucial in spectral clustering, similarly holds for them.
In addition to its significance as a theoretical analysis of the t-SNE algorithm a popular method for dimen- sionality reduction and visualization-tie there is practical utility that the ODE we derive has a closed-form solution. This means that by performing eigenvalue decomposition of the adjacency matrix in advance, we have the advantage of obtaining embedding results at any desired point in time.
13

---

Published in Transactions on Machine Learning Research (06/2025)
Ik=27, 1=2700)
(=1400)
1=280.0
ARR/ / Trustworthiness for MNIST dataset
X00
ARR frustworthiness
Method & Metric
GD ARR
330
GD TW
MM ARR
<VC MV
-30 KM 00 100 05
28 4001 Hi on 10
MM TW
GD (c=200)
MM (K=200)
MAG (x=200)
NAG ARR
NAG TW
2390
33/8
-3-
-34
Figure 2: (Jeft): transition of ARR and Trustworthiness (TW). (right); clustering result for MNIST dataset with ARR=0.01 using three ODEs GD, MM, and NAG (right-top), and result after the following embedding stage (right-bottom)
While this framework currently relies full eigendeccompssiioon of the graph Laplacian, which incurs O(n²) complexity, we emphasize that this step is used primarily for thcoretical analysis and computation of diagnostic metrics such as ARR. For large-scale applications, we anticipate that approximate met such as the Lanczos algorithm or randomized eigendecomposition could reduce the computational burden to around O(n²), especially for sparse graphs. Exploring these approximations would a promising direction for extending our framework high-dimensional datasets at scale.
We discuss the approximated formula for t-SNE only concerning the EE stage but do not address the embedding stage. Theoretical analysis of the embedding stage will complement work and further deepen the understanding of t-SNE. Additionally, we incorporate bandwidth into the perplexity but only examined specific case. broader analysis of different scenarios and their implications on perplexity would provide a more comprehensive understanding.
UMAP (Ghojogh et al., 2021 has been frequently used alongside t-SNE. There are studies discussing the relationship between t-SNE and UMAP (Damrich et al. 2023 ; hence, analyzing UMAP through dynamical systems-based approach is intriguing direction for future research.
Acknowledgement
Part of this work is supported by JSPS KAKENHI (23K24909 and 25H014944). Finally, we express our special thanks to the editor and anonymous reviewers whose valuable comments helped to improve the manuscript.
References
Sanjeev Arora, Wei Hu, and Pravesh K. Kothari. An analysis the t-SNE algorithm for data visualization. In Sébastien Bubeck, Vianncy Perchet, and Philippe Rigollet (eds.), Proccedings of the 31st Conference On Learning Theory, volume 75 of Proceedings of Machine Learning Research, pp. 1455-1462. PMLR, 06 09 Jul 2018. URL httpps://rroceedings. mlr wwtp:pppss///sss.aaaalaao
http://cccccaagngn.pprrrreeeaaaa htmml
14

---

Published in Transactions on Machine Learning Research (06/2025)
Antonio Auffinger and Daniel Fletcher. Equilibrium distributions for t-distributed stochastic neighbour embedding. arXiv preprint arXiv:2304.03727, 2023.
Sivaraman Balakrishnan, Min Xu, Akshay Krishnamurthy, and Aarti Singh. Noise thresholds for spec- tral clustering, In John Shawe-Taylon Richard S. Zemel, Peter L. Bartlett, Fernando C. N. Pereira, and Kilian Q. Weinberger (eds.), Advances in Neural Information Processing Systems 24: 25th Annual Conference on Neural Information Processing Systems 2011. Proceedings of a meeting held 12- 14 December 2011, Granada, Spain, pp. 954-962, 2011. URL httppp://prooeddings...rr.spssps wttppppppsccpppper//.eeea.ap httppp/pp////hhhhhh
html
Mikhail Belkin and Partha Niyogi. Laplacian eigenmaps and spectral techniques for embed- ding and clustering. In Thomas G. Dietterich, Suzanna Becker, and Zoubin Ghahramani (eds.), Advances in Neural Information Processing Systems 14 [Neural Information Processing Natural and Synthotic, NIPS 2001, December 3-8, 2001, Vancouver, British Columbia, Canada]: : pp. 585 591. MIT Press, 2001. URL htttp::/prroceedings www..nuuii...om wttp://pp...ccpppp//////////aaaa
Systems:
11067/0900c11300cca00ccccoeeee Abstrccctbbbrtttt.... html
Mikhail Belkin and Partha Niyogi. Laplacian eigenmaps for dimensionality reduction and data representation. Neural Comput., 15(6):1373 1396, 2003. doi: 10.1 10.1162/089996660332277777 URL https://doi. org/10.
1162/0899660333788177
Jan Niklas Böhm, Philipp Berens, and Dmitry Kobak. Unsupervised visualization of image datasets us- ing contrastive learning. In The Eleventh International Conference on Learning Representations, ICLR 2023, Kigali, Rwanda, May 1-5, 2023. OpenReview.net, 2023. URL sttp:://opemneeviee..eee//tt
n 12Hmmtt
T. Tony Cai and Rong Ma. Theoretical foundations of I-SNE for visualizing high-dimmensionall clustered data. J. Mach. Learn. Res., 23:301:1 13:101:1-301:54, 2022. URL http://jmlr http://iwi..oorgppppesssss////// orgg/pppers//22221-054444 html
Miguel Á. Carrcira-Porpinánn The clastic embedding algorithm for dimensionality reduction. In Johannes Fürnkranz and Thorsten Joachims (eds.). Proceedings of the 27th International Conference on Machine Learning (TCMT,10), June 21-24, 2010, Haifa, Tsrael, pp. 167-174. Omnipress, 2010. URL https://imml.
cc//oofference////////peerr////e 11..ppeerrer1.............
Pin-Ya Yu Chen, Baichuan Zhang, Mohammad Al Hasan, and Alfred O. Hero III. Incremental method for spectral clustering of increasing orders. CoRR, abs/1512.07349) 2015. URL http://arxiv. org/abs/
1512.07349
Sebastian Damrich, Jan Niklas Böhm, Fred A. Hamprecht, and Dmitry Kobak. From t-SNE to UMAP with contrastive learning. In The Eleventh International Conference on Learning Represcntations, ICLR 2023, Kigali, Rwanda, May 1-5, 2023. OpenReview.nel, 2023. URL ftttp:://openreview.eeet www/net ?id=
B8a1FcYOvi
Aaron Defazio. On the curved geometry of accelerated optimization. In Hanna M. Wallach, Hugo Larochelle, Alina Beygelzimer, Florence d' Alché-Buc, Emily B. Fox, and Roman Garnett (eds.), Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8- 14, 2019, Vancouver, BC, Canada, pp. 1764 1773, 2019. URL. httpp::/ wtwp:////33000-t-aarrrt.... html
proceedings.o
Takanori Fujiwara, Jia-Kai Chou, Shilpika, Panpan Xu, Liu Ren, and Kwan-Liu Ma. An incremental dimen- sionality reduction method for visualizing streaming multidimensionall data. IEEE Trans. Vis. Compul. Graph., (6(1))418-42)) 428, 2020. doi: 10.1199/TVCG..219.2200 2934433. URL https://doii org/10. wttgg//...///////........ 1109/TVCG. 2019. 2934433
Benyamin Ghojogh, Ali Ghodsi, Fakhri Karray, and Mark Crowley. Uniform manifold approximation and projection (UMAP) and its variants: Tutorial and survey. CoRR, abs/2109.02508, 2021. URL https:
15

---

Published in Transactions on Machine Learning Research (06/2025)
Shin-itiro Goto and Hidcitsu Hino. Fast symplectic integrator for nesterov-type acceleration method. Japan Journal of Industrial and Applied Mathematics, 42(1))331132) -32, 2025. URL
813160-024-00880-4
Yunhui Guo, Haoran Guo, and Stella X. Yu. CO-SNE: dimensionality reduction and visualization for hyperbolic data. In IEEE/ /CF Conference on Computer Vision and Pattern Recognition, CVPR 2022, New Orleans, LA, USA, June 18- 24, 2022, pp. 11-20. IEEE, 2022. doi: 10. LLLLos/CVPR2268.222222.2.
URL https://dol. http://l....rrr....1/////////e// 1119///P552688..... 2022.0002
Zhiyong Hao, Yixuan Jiang, Huihua Yu, and Hsiao-Dong Chiang. Adaptive learning rate and momentum for training deep neural networks. CoRR, abs/2106.11544 2021. URL
11548
Lulu He, Jimin Ye, and Jianwei E. Nonconvex optimization with inertial proximal stochastic variance reduction gradient. Inf. Sci., 648:119546, 2023. doi: 10.1016/..INS.2233.11000 URL https://dol: org/
10. 10.1016/ 2023.202..00 119546
Gooffrey E Hinton and Sam Roweis. Stochastic neighbor embedding. In S. Bocker, S. Thrun, and K. Obermayer (eds.), Advances in Neural Information Processing Systems, vol- ume 15. MIT Press, 2002. URL https:// http://roooddngnn.urrr.p....... http://appppeeleeeeppeeeeeeeeeee fit.ss//ppper222//Ei/2 pdf
Lawrence Hubert and Phipps Arabie. Comparing partitions. Journal of classification, 2:193 218, 1985.
Seonghyeon Jeong and Hau-Tieng Wu. Convergence analysis of t-sne a gradient flow for point cloud on manifold, 2024.
L.. Jolliffe. Principal Component Analysis. Springer Verlag, 1986.
Samuel Kaski, Janne Nikkilä, Merja Oja, Jarkko Venna, Petri Törönen, and Eero Castrén. Trustworthiness and metrics in visualizing similarity of gene expression. BMC Bioinform., 4:48, 2003. doi: 10.1186/ 1471-2105-4-48. URL
Dmitry Kobak and George C Linderman. Initialization is critical for preserving global data structure in both t-sne and umap. Nature biotechnology, 39(2):156-117, 2021.
Dmitry Kobak, George C. Linderman, Stefan Steinerberger, Yuval Kluger, and Philipp Berens. Heavy- tailed kernels reveal a finer cluster structure in t-SNE visualisations. In Ulf Brefeld, Élisa Fromont, An- dreas Hotho, Arno J. Knobbe, Marlocs H. Maathuis, and Céline Robardet (eds.), Machine Learning and Knowledge Discovery in Databases - European Conference, ECML, PKDD 2019, Würzburg, Germany, September 16-20, 2019, Proceedings, Part 1, volume 11906 of Lecture Notes in Computer Science, pp. 124 139. Springer, 2019. doi: 978 46150-888
8. URL https://doi http:///ooooorg/..1101 org/10.1007/
Nikola B. Kovachki and Andrew M. Stuart. Continuous time analysis of momentum methods. J. Mach. Learn. Res., 22:17:1-17:40, 2021. URL http:://malr..rrr/ppprs//////
html
Jesse H. Krijthe. Rtsne: t-Distributed Stochastic Neighbor Embedding using Barnes-Hut Implementation, 2015. URL http://://ggttuu.b..o////rrrr//e
R package version 0.17.
Solomon Kullback and Richard A Leibler. On information and sufficiency. The annals of mathematical statistics, 22(1):79-50) -87999966111 1951.
Pierre Lambert, John A Lee, Edouard Couplet, and Cyril De Bodt. Nesterov momentum and gradient normalization to improve t-SNE convergence and neighborhood preservation, without carly exaggeration. ESANN proceedings, 1:1, 2023.
George C. Linderman and Stefan Steinerberger. Clustering with t-sne, provably. SIAM J. Math. Data Sci., 1(2):313 332, 2019. doi: 10.1137 2.1137/18M1216133 URI. httpp:://aao..rrr//1111/////0010
16

---

Published in Transactions on Machine Learning Research (06/2025)
George C. Linderman and Stefan Steinerberger. Dimensionality reduction via dynamical systems: The case of t-sne. STAM Rev., 64(1):153-17) 157118 2022. doi: 10.1137/21M1446666, URL /ttp///6....g//11111/////
21m 1lm14465 446769
Leland McInnes and John Healy. UMAP: uniform manifold approximation and projection for dimension 2018. URL
reduction. CoRR,
1102//211220006666 03426
Kevin R Moon, David van Dijk, Zheng Wang, Scott Gigante, Daniel B Burkhardt, William S Chen, Kristina Yim, Antonia van den Elzen, Matthew J Hirn, Ronald R Coifman, Natalia B Ivanova, Guy Wolf, and Smita Krishnaswamy. Visualizing structure and transitions in high-dimensional biologi- cal data. Nature biotechnology, 97(12):1482-144 December 2019. ISSN 1187-0156,1546-111 1696. doi:
00.1088//441587--19-03366
Ryan Murray and Adam Pickarski. Large data limits and scaling laws for tsne. CoRR, abs/ 2212221111...313 2024. doi: 10.48550/ARA ARXIV 1310. URL httts://doi.org:: 00tg///..8550//800 w...485555Xii11 /arviv.Xiivv///// 2410 13063
Yurii Nesterov. A method for unconstrained convex minimization problem with the rate of convergence oL1/k²). Soviet Mathematics Doklady, 27:372 376, 1983.
Frank W. J. Olver and Leonard C. Maximon. Bessel functions. In Frank W. J. Olver, Daniel W. Lozier, Ronald F. Boisvert, and Charles W. Clark (eds.), NIST Handbook of Mathematical Functions, pp. 215 286. Cambridge University Press, 2010.
Antonio Orvieto and Aurelien Lucchi. Continuous-tine models for stochastic optimization algorithms. In Advances in Neural Information Processing Systoms 32, pp. 12610-12622. Curran Associates, Inc., 2019.
Nicola Pezzotti, Boudewijn P. F. Lelieveldt, Laurens van der Maaten, Thomas Höllt, Elmar Eisemann, and Anna Vilanova. Approximated and user steerable t-SNE for progressive visual analytics. IEEE Trans. Vis. Comput. Graph., 23(7):1739-13 1752, 2017. doi: 10.1109/TVCG, 1LL10/TVCC..011.55775555 URL https: 1109/TVCC. 2016. 2570755
Zana Rashidi, Kasra Ahmadi K. A., Aijum An, and Xiaogang Wang. Adaptive momentum coefficient for neural network optimization. In Frank Hutter, Kristian Kersting, Jefrey Lijffijt, and Tsabel Valera. (eds.), Machine Learning and Knowledge Discovery in Databases European Conference, ECML PKDD 2020, Ghent, Belgium, September 14-18, 2020, Proceedings, Part II, volume 12458 of Lecture Notes in Compuler Science, pp. 35 51. Springer, 2020. doi: 10. 110007078---0001 3-030-67661-22- 3. URL attps://doi. o0./10.
Sam T. Roweis and Lawrence K. Saul. Nonlinear dimensionality reduction by locally linear embedding. Science, 290(5500):2323 2326, 2000. doi: 10.1126/ (0.1126/science. 2tt.ce...0.550.2222222. URL
org/doi/abs/10..ooo 11tp:.......1aaeeece.e 290. 5500. 2323
Martin Skrodzki, Nicolas F. Chaves-de-Plaza, Klaus Hildebrandt, Thomas Höllt, and Elmar Eisemann. Tuning the perplexity for and computing sampling-based based t-SNE embeddings. CoRR, hts//2008.1500 2023. doi: 10.48550/ ARXT 200...23..111 15513. URL ltttts::///i..ooo //tp::///oooooo/////// wtw./////...///.///////// 08550/aaXix. 2308 15513
Martin Skrodzki, Hunter van Geffen, Nicolas F. Chaves-de-Plaza, Thomas Höllt, Elmar Eisemann, and Klaus Hildebrandt. Accelerating hyperbolic t-sme. CoRR, abs/2401 13700 13708, 2024. doi: 10.48550/ A 2401.
13708. URL https:://oi. 10. 4w0544885500i.../ 13708
Fan Wei Lee Wenke Prodromidis Andreas Stolfo, Salvatore and Philip Chan. KDD Cup 1999 Data. UCI Machine Learning Repository, 1999. DOI:
Salvatore J. Stolfo, Wei Fan, Wenke Lee, Andreas L. Prodromidis, and Philip Chan. KDD cup 1999 data. https://doi org/10.
10. 24432/C51C7 December 1998. URL https://doi. / ottp://0....rr///..0.4433//11111 Ac- cessed on 2024-12-05.
24432/C51C7N
Weijic Su, Stephen P. Boyd, and Emmanuel J. Candès. A differential equation for modeling nesterov's accelerated gradient method: Theory and insights. J. Mach. Learn. Res., 17:153:1-153:43, 2016. URL
http://mmr..orr//ppprss///////// 1ttmmmmmm....
17

---

Published in Transactions on Machine Learning Research (06/2025)
Ily Sutskever, James Martens, George E. Dahl, and Geoffrey E. Hinton. On the importance of initializa- tion and momentum in deep learning. In Proceedings of the 30th International Conference on Machine Learning, ICML 2013, Atlanta, GA, USA, 16-21 June 2013, volume 28 of JMLR Workshop and Conference Proceodings, pp. 1139-1147. JMgR.org, 2013. URL htLp:tt::// htttp::ppprocediinggg/ mlr www...reessvv.... preess/v28/utskeeerr13 sutskever13. html
P. Szymański and T. Kajdanowicz. scikit-based Python environment for performing multi-label classifi- cation. ArXiv e-prints, February 2017.
Mahbod Tavallaee, Ebrahim Bagheri, Wei Lu, and Ali A. Ghorbani. A detailed analysis of the KDD CUP 99 data set. In 2009 IEEE Symposium on Computational Intelligence for Security and Defense Applications, CISDA 2009, Ottawa, Canada, July 8-10, 2009, pp. 1-6. IEEE, 2009. doi: 10.1109/CISDAA /OO/CCCAA...9..5555535 URL https: //op:a://doi. org/10. 1109/CISDA, 2009. 5356528
J.B. Tenenbaum, V. de Silva, and J.C. Langford. A global geometric framework for nonlinear dimensionality reduction. Science, 290(5500):2319-2555 23233, December 2000.
Warren S Torgerson. Multidimensional scaling: 1. theory and method. Psychometrika, 17(4):401 419, 1952.
Laurens van der Maaten. Barnes-hut-sne. semanticscholar wwwp:oooorot.t
CoRR, abs/ 130/133 2013.
URL https://api.
senaatttccccblaa..rr///oppoposss 2089155226
Laurens van der Maaten. Accelerating t-sme using tree-based algorithms. J. Mach. Learn. Res., 15:3221-3245 2014. URL htt.tpp://app..eeaattccccoooooo http://pwppaapparooooooooooorror :1: 14618636
org/CorpusID:
Laurens van der Maaten and Geoffrey Hinton. Visualizing data using t-SNE. Journal of machine learning research, 9(11), 2008.
Jarkko Venna and Samuel Kaski. Neighborhood preservation in nonlinear projection methods: An exper- imental study. In Georg Dorffner, Horst Bischof, and Kurt Hornik (ods.), Artificial Neural Networks ICANN 2001, International Conference Vienna, Austria, August 21-25, 2001 Proceedings, volume 2130 of Lecture Notes in Computer Science, pp. https://aoi.org// 0ttp://00..000000-004444----0000
Springer, 2001. doi: 10. 44640004468-0,, 68. URL
Max Vladymyrov and Miguel Á. Carreira-Perriiñánn Fast training of nonlinear embedding algorithms. ArXiv,
abs/ Abs/1206.4646 2012. URL https://ppp..smmaatiiccooo...e http://ccchcooaaoooorrrorppsspee Ortg/oorpsID::eee.cooo
Max Vladymyrov and Miguel Á. Carreirra-Peppiãánn Linear-time training of nonlinear low-dimensional embeddings. In International Conference on Artificial Intelligence and Statistics, 2014. URL https:
//api.semantisccho.... org/CorpusID:1083216
Uhike von Luxburg. A tutorial on spectral clustering. Stat. Comput., 17(4):395-411 2007. doi: 10.1007/ $11222-007-9033-2 Z. URL htte:://0.....////00000///112210 -033--
George Neville Watson. A treatise on the theory of Bessel functions, volume 2. The University Press, 1922.
Andre Wibisono, Ashia C. Wilson, and Michael I. Jordan. A variational perspective on accelerated methods in optimization. Proceedings of the National Academy of Sciences, 11355 2016. ISSN 0027-8424. doi: 10.1073/pnas, 1..1077/pmaa..61444111133 1614734113.
Ashia C. Wilson, Benjamin Recht, and Michael I. Jordan. A lyapunov analysis of accolerated mothods in optimization. J. Mach. Learn. Res., 22:113:1-113:3, 2021. URL sor.tccco...mmmtiiccccccco...... oorgg/
CorpusID:235341488
Zhirong Yang, Jaakko Poltonen, and Samuel Kaski. Scalable optimization of neighbor embedding for visual- ization. In International Conference on Machine Learning. 2013. URL semanticscholar.
14514454
18

---

Published in Transactions on Machine Learning Research (06/2025)
Ryeongkyung Yoon and Braxton Osting. A dynamical systems based framework for dimension reduction. CoRR, 08155,, 2022. doi: 10.48550/ ARXTV22040000015555 arxiv. 08155
URL. https:://db..org/10.448000
Yi Yu, Tengyao Wang, and Richard J Samworth. A useful variant of the davis-kahan theorem for statisticians. Biometrika, 102(2):315-323, 2015.
Jingzhao Zhang, Aryan Mokhtari, Suvrit Sra, and Ali Jadbabaie. Direct runge-kutta discretization achieves acceleration. In Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Biande Cesa-Bianchi, and R. Garnett (eds.), Advances in Neural Information Processing Systems 31, pp. 3900 3909. Curran Associates, Inc., 2018.
Yuansheng Zhou and Tatyana Sharpee. Hyperbolic geometry of gene expression. Iscience, 24(3), 2021.
Ye Zhu and Kai Ming Ting. Improving the effectiveness and elficiency of stochastic neighbour embedding with isolation kernel (extended abstract). In Lud De Raedt (ed.), Proceedings of the Thirty- International Joint Conference on Artificial Intelligence, IJCAI-22, pp. 5792 5796. International Joint Conferences on Artificial Intelligence Organization, 2022. doi: 10.24963/ 2022 2022 /812. URL https: //doi. org/10.
24963/1jcai. 2022/812 Journal Track.
19

---

Published in Transactions on Machine Learning Research (06/2025)
A Derivation of the t-SNE gradient coefficient
For the sake brevily, we omit the iteration index k; we use ye, Piys Oijs instead of yi Pij and lij" The contents basically follows Appendix A in van der Maaten & Hinton (2008).
(k) (2)
(k)
t-SNE minimizes the KL divergence between the joint probabilities Pa in the high-dimensional space R and the joint probabilities 4ij in the low-dinensional space Rd: They are defined in Eq. (1,, P.j (ptj + Pyte)/ 2n and Eq. (2, respectively. Cost function C' is defined with KL divergence Kullback & Leibler (195]
C (LLII))
Σ Pij log
Pij hi
Σ
(pgy log Pij - Pig log 9:))
(A.1)
For the purpose of the derivation briefly, we set two variables dip and Z:
di = 11 Myll
(A.2)
+ dz
(A.3)
Then, with the chain rule, the gradient of the cost function with respect to yi is given by
ac
8C
0C
Σ
+
(yi w;)/d:;
(A.4)
Oy:
Od,j
dd.z
aC
= 2 Σ
(ys U;)/
(A.5)
dd.j
Pij are defined with high-dimensional data points, and they're irrelevant to ys or dig that is 0d.5
Opac
0 for
Yk, I E [n]. Then, we have
OC
i llg qui
Pkt
(A.6)
Od.j
Od.j
A
alog Oki 7 log Z)
Pit
(A.7)
Öd,j
@(11 + 11)
1 OZ
Σ pst
(A.8)
que Z
Odg
Z ddis
kyf
The gradient term a((++d,)) 1 Odiy
is only nonzero when = ind j. Then,
ac
Pij (1 qu
d 2d
Σ keft
Pkl
(1 + de, -
(A.9)
0d,j
Z
Note that ExSt Pii 1 and we have
ac
21:s(1 -/ 21, (1 + dl.) -ldj 2(pis -2)) d dij.
(A.10)
0dg
Finally, we have
ac Oy: 4 Σ
Pij Øs
1+ 11: y, 12 (v: 4z).
(A.11)
In this formula, we aim to emphasize the effect of the original data distribution using Pij. We formally substitute Pty with apy and define the term 1t apg -4
11 1 as (a)) This leads to the expression in Eq. (4).
20

---

Published in Transactions on Machine Learning Research (06/2025)
B The effect of accelerator parameter
In this section, we examine the effect of the acceleration parameler a separately for GD and MM/NAG.
B.1 The case of GD
First, consider the update equation (B12) for GD. The low-dimensional representations y and yj (k) (i + j)
show that [k+1) is calculated using y/ (k) and gradient term. The sign of opij -9g (k) inffuence. We define aptj the attractive force and lij as the repulsive force. When apš - 9ij > 0,
4g determines its directional
mcaning the attractive force is stronger, y/ (k+1) moves closer to Vj (k) compared to y. Conversely, when (k) apij Pig < O, the repulsive force dominates, causing Vi (k+1) to move further away from v,th Figure 3 illustrates this siluation.
(k)
Vi
(k+1)
Vi
(6) h Σ
Opsj
Dij
(1)
(k)
Di
(1)
(B.12)
1 1101
(k)
(1c
ISiJSn, sfy
y 2
yj (10) (K)
Yi
(k) (repulsive force)
in case that (attractee force) qij
In case that (atttttrace force)
< <q (repulsive force)
(k)
00
(k+1)
V/ (k)
Figure 3: Diagram of the update equation for y/
(k)
Therefore, the larger the value of 1, the stronger the attractive force among nodes becomes. As result, this generally facilitates the formation of clusters.
B.2
The case of MM and NAG
As the update equations (10 for MM and (19) or 20 for NAG show, yi (k+1) calculated using both the gradient term and the momentum term. Accelcration achieved by incorporating not only the current gradient but also past updates. Even when using MM or NAG, the amplification of the attractive force with tends to occur.
21

---

Published in Transactions on Machine Learning Research (06/2025)
In summary, when 1, the attractive force generally becomes stronger, making it casier to form clusters across all optimization methods GD, MM, and NAG. This suggests that larger can accelerate the clustering process.
C Asymptotic behavior of the update equation
C.1 Proof of Theorem
This theorem Cai
Ma 2022) based without. momentum term in 43) and can be extended even
with momentum term, because the proof docsn't refer the update equation (5. but just the definition of
$ (a) in Eq. 4
C.2 Proof of Proposition 4.2
This proof proceeds in the same manner as the proof of Proposition 3 Cai & Ma 2022 (k) || hs! 1 lllll 1lo0 + mly+
Note that
(k+1)
(k)
yt
(k-1)1 11 for any k 2 1, where
n h
||T As I My 110
Σ
S (a:
+
|
11 2h Σ
S (a)l-
(k)
j=1
ft
For the second term, we have
h Σ 1 (a)l < 1 s nn(aa||Pll> + [[]])a)) S nha Pllx + hs1 + 1) ()) n
where the last inequality follows from Eq. (40) in Cai & Ma 2022 (k)
Intermediate variable Z defined
Ents( (1 + llu (k) ") 12 and we can deduce it with | (k) 10 S 1/Z, Z 2 n(n 1)/(1 + v)()) and
Q 10 $ (1 + 17 (z)) /n(n 1). Then we have
h(1 + ()))
1k+1)
(1 2nha Pl 00 +
11x (1) m m2+11)) (8+1) (1y 1100 m (x+ 11x (1))((/-)11 1100
Vi
11
h(1 + 17 (k)
(1 + 2nho P| +
m(++1)) 11///01100 + m (x-1)) |11 (k- 11- 11000
n1
1(1+ + (k)
1 + 2nha||
+
(++1) (11y م| + ! (1-1) |1اo)
2h(1 + 17(0)
2 + 4nha P| +
2m(
max max (*)) وحاا lly (1)) bا 1e
n 1
2h(1 + (T))
( 1|x0 < (2 Annha||/llllle
2m (+t1))
max (1ly 1mc 1y((~1111>>>
1E(2
Introducing Tn = nhai P Pooloo + A we have
( 1oo0 (2 2m (k+1) + Cra)
max
(1y
11,0, 1ly/ 11)(--1) [10)
(C.13)
[E[3]
and for any k 2 2,
(x 11000 S (2 2m (k+1) + Cra)
max f (1) 1loo, |1y/ (0) 1100 lEp2
(C.14)
Lemma C.I. For Ak 2 2, (k) and max (1ly) (k)
loc: Ily!
( bock bounded by a constant.
22

---

Published in Transactions on Machine Learning Research (06/2025)
Proof. In case of k= = 2, Eq. (C.14) shows that
ly 1oc (2 (m( Crn) max (1ly)" 1boo, 11y/"1100)
(C.15)
The condition (TI) states that max1e121 |1y 1o, (lyl1 1lo) is bounded. According to C.14).
( - 4(2 + (m) + Cr,)²
max 10 11y) 12
(C.16)
[E[2]
is also bounded. Then, let's say that the statement holds for k, i.e. both Iyy (k) vm and [ly[++D) bl are bounded. With incquality CC.13). Ily/
(xy(++))) lloe is also bounded. Assuming on (1) by condition (T1)
v (++1) < 4 max
1 < max {lly (k+1) 1E(2
1 V2
(k+1)
1
(C.17)
1E(M)
< 4(2 + 2m (++1) + Cra)² max
(Ily
( (k-1) 12at O(1),
1 I1y
(C.18)
[E(2]
we have (8++1) is bounded. By induction, the statement holds for Ak: 2 2.
As long as k(n) with kv, O(1) by condition (TI), we have
Ily ( /l///
max (Ily) (2) Ioo, Ily lyy" =
(C.19)
[E
or,
diam {,y (1:)
}oe(m)
maxie[n) 13/1 (A)
O(1).
(C.20)
max (lly]
(1)
bl lly (0)
bloc
max (lly|" bo 11y (0)
11
It proves the statement.
D Derivation of ODEs
We proceed to concretely examine the following two cases: MM and NAG.
D.1 Proof of Lemma 5.1
The assumption
(khh)) is justfied with the following salatement, which is shown similar to Proposition 8 in Cai & Ma (0022). We consider the momentum terms to be constant mn( 1), i.e m m(+1)) for vk. So, we can rewrite the formula (8) as follows:
y (++1)
xe II - hL(PP H,y + m(y - yi ((=-1)),, L [2].
(D.21)
For 1 E [2], let [8] (k)
1220 be the sequence defined the iteration with the update formula:
((++1) [I, hL(oP H,))y((*) + m()(") s((-1)), k21
(D.22)
with initial values (1)
y()) y = y Introduce the assumption Y((i) Rs (k) for some smooth curve
(0)
(tt defined for t 2 0. Set t kh. Then, the update formula (D.22) reduces
Y,(t + h) Y(()) - ML(aP H,) + ml (Y,()) Y((l hh))
(D.23)
By dividing by > O, we have
Yi(l + h) Y(()) -L(aP HY))((1) + Y((I) Y,(L - -h)
(D.24)
h
h
From the definition of the derivative & YY limm (((1+))-Yii -im (( have
15
23

---

Published in Transactions on Machine Learning Research (06/2025)
D.2 Proof of Proposition 5.2
For any Ago, we set e = y/ (k(h),, L= L(oP, - -H,) for brevity.
Lemma D.1. Let (t( the formula in Proposition 3.1 and a symmectic matric L has an eigendecomposi- tion,
L L=
(D.25)
where A diag(Ai.... An), A S An and Q E O(n). For t >>
(i(())22 - exp (- 11-m m-m/ ) |Y(())22
ty1
(D.26)
Specifically, 0 < < 1 and the condition (I1) holds, [f |l is bounded. Proof.
With the eigendeccomposition L QAQ-! 1, we have
)| |Y(())222
(t)()) <
exp
1- m
(D.27)
exp
t
A Q- (
(D.28)
1 1- m
exp
A
(D.29)
1-m m
S exp
(-1--,)
[Y)l2)
(D.30)
m
Using Taylor expansion for some & E (kh, (I: + 1)h), and Lemma 5.1
we have
Yi((k + 1)h) = Yi(kh) + (kh)
122 d Y
(D.31)
di
2
d12
h
h22
1
Yi(kh) 1- m (kk)
Y(S)(
(D.32)
1 m
We set those formulac into
ek+1 = Yi((k + - -5())s(ttt)))
(D.33)
ck+1 (Yi(kh)
h I-m m (kh) + [/2 h
2
)())
(D.34)
2 1 ml
- m)) - -ny!! - -my!* (1-)
(D.35)
With (kh) = (A) + ck: we have
Cx+1 Cx m
hm
(kh) + m(s((t1)) -s(th) + 2 1 m
1 h
2
L2Y,
(D.36)
Due to triangular inequality, we have
|et+illz ll2 S ||ex l12 +
hm
| 122 mlls((-1)) -)(²² 1|2 +
(k-1)
1 1 h 1 (
(D.37)
1 m
2
Summing up on the index, we obtain
|ex+1ll2 < ||eoll2
Σ 1=0
hm |||||! 2l2 + ml]s((-11)
1|2
( h
(1_mm)² " (ILI))X)((())
(D.38)
2
24

---

Published in Transactions on Machine Learning Research (06/2025)
Thanks to initial condition in Lemma 5.1 en -5( Y:(0) (0, ;)", Proposition 4.2 says that have constant C1 > 0 such that
max; (3)
( Ci.
(D.39)
Additionally, with respect to the term || (1) ll2, consider the update equation (10
8 - 51 ((1-1) m(F
1 y )(1 - ALA! 1
(D.40)
We have
1 - 1|2 S mlly 5 + hC1 h
(D.41)
By recursively applying these incqualitics, we obtain
| ()1)) (1)) 1|2 S (m' + m +... + m + 1)/Ci|LL ||)12))
(D.42)
1 - m²
hC1
L (Y())
|L| (0)l22
(D.43)
m
With Lemmaa
D.1 and (11),
Y (1)
Eo
(0)l2 < exp (-1- m
1-
S
exp
(D.44)
m
We denote the right-hand side as C2 > 0. Under these conditions, we have
lek+ l|2 (o)l2
2hm + 1)C1
++1 1 h
2
v
L +
|2
(D.45)
1 1-m
2
-
This proves the statement.
D.3 Proof of Lemma 5K3
Since
L(P oH, L(aP H,) 1
1
In n(n 1)
1 1
(D.46)
the eigenvect.ors L(aP H.) and L(aP) share the eigenvectors. L(P aH,)) = 1 u.
Also, if we decompose eigenvector as
L(aP H,)u, = d.44.
(D.47)
We have 01 adi,o; oj ad T for 2 2 n, which corresponds to 14
D.4 Proof of Lemma 5.5
Remember that both (L(()))) S (41(L(P))) i [n- 1] and AR+1 (L(()) is the first positive eigenvalue, )R(L(P")) = 0, A (L(")) 0.
i.e.
(D.48)
With Weyl's inequality
(L((P)) - (L(L)P)) - [(E)]l:
(D.49)
where PPPPP- P When R in incquality. D.49 we have
(L [(E]]]
(D.50)
The condition (T2.M) says that - L(B) o(1) and it's equivalent that
We N,3N E s.t. n n N t m
[(E)] <
(D.51)
25

---

Published in Transactions on Machine Learning Research (06/2025)
If we set (n-1)(1-m) have
1
adR(L(P)) n-1 a
[BB]
(D.52)
m)
(D.53)
t
n -
x
m)
0.
(D.54)
(n 1)(1 - m)
n
By contrast, when i= R + 1 in Ineq. (D.49
A AR+1 ( |L(E)]]
(D.55)
So, we have
A 2 (L+1)((P))) - |[(EE]|-
(D.56)
If we set = e= (1- ----) x
I-m and the condition (T2.M) such as AR- (L=((((P)) > (la)-1
a)n+i(L(P)) n
2 I az L(E)]
(D.57)
n
(1 m)c
1
(L
(D.58)
T-
m x (1
> a
x
= 0.
(D.59)
to
1
1-m m
n n-1
Note that the condition (TIM) t = o(n) assures that = (1 n-1 x 1-m
> for n>
and 0 < m < 1.
E
Well-conditioned matrix
As discussed in section 5.1 we're interested the behavior of (t for > 1. The sign of the terms aλ; n-I decides it; if negative, corresponding components are amplified: if positive, they are suppressed as noted in 17
This scenario, where aλ; - < 0 for 2 i R, indicates that the i-th eigenvalue X is close to 0 under a > 0. Although the matrix L(P), by definition, has one connected component, suppose there exists another matrix P sufficiently close to P such that its Laplacian L(P") has R connected components, implying the dimension its mull space is R. As discussed in CCai & Ma 2022). it natural to pick up Laplacian mull space. We call the adjacency matrix p as "well-conditioned" connected components.
its associated weighted graph has 2 2
Proposition E.1. (Proposilion 6 in (Cai & Ma 2022) Let A E Rax* be symmelric and well-conditioned matria. Then, the smallest eigenvalue of the Laplacian L(A) is 0 and has multiplicity R, and its associated eigenspace is spanned by {0....., OR), where for each r E [R]
1/ 1/4 if node belongs component i,
00
(E.60)
0, otherwise,
and ny is the number of nodes related to y-th connected component. Also, the rull space of L(A) is spanned with these basis
0
n
Ang
:
(F.61)
n22
/RR
0
1
where In, is length n, column vector with ones, and 0 is vector with appropriate length.
When the data (Xiheen are well clustered and bandwidths T; are appropriately selected, in (Balakrishnan et al. 2011). it is shown that we have good approximated matrix of P. Such well-conditioned matrix P is obtained concretely in section g Such matrix P is theoretically used in implicit regularization theorem 6.5
26

---

Published in Transactions on Machine Learning Research (06/2025)
F
Derivation of ODE for NAG
F.1 Proof of Lemma 6.1
With two equations
19 and
20 we derive another ODE. Applying a rescaling, we have
(k++1) 1)
(R) Vh
k:
(K)
y y(l-1)
ViL(aP - -H,) H
(F.62)
h: + Vh
Introduce the [assumption] 1/1h) Y(kVh) for some smooth curve Yi(t) for O. Set K t/vh. Then,
as the step size h goes zero, Yi(t) = y ((1/vh) spth and Yi(t + vh) s (((t/B)//)) y1 (k++)) Taylor expansion gives
(y( y D/V + 2 + (y(t)) yi ( - -
and VhL(aP -H, - VhL(aP ( oVh)). So, F.62 is written as Y((0)+ 2 v = (1-3 3v (i((t) 1 2 +(00)++ - vTL((PP
- H (.(t) +o(Vh).
By comparing the coefficients of Vh, we achieve
Y,(() + - 1
+ (()))) =0 E [2].
(F.63)
Note that the first initial condition is Y,(0) = y Also, taking k in F.62) we have
(x
//T/ -y VTP (,, o(1)
yi
Then, the second initial condition just Y,(0) (0, (O.....0)) ER.
F.2 Proof of Proposition 6.2
If we define a function f: Rd R. as f(y) by L(aP- H.Jy for y R", the function f differentiable and fz)) L(aP H,yy. So, we have
[V)(2) V f(y)l2 - on 0, y112.
(F.64)
where nn is the maximum cigenvalue of L(oP H.), especially on -ipsccitz function. So, apply Proposition 2 in Su el al. 2016), which shows our stalement.
F.3 Proof of Lemma 6.3
Similar to Lemma 5.5 we deduce the results since have 1 L(P P*)]| = o(1), that is
VE > 0, EN E N s.t. n n N +
<e,
(F.65)
(L+1(((P))) > (Pa)-1 which are both (T2.N) and o(n as in (TI.N).
F.4 Proof of Proposition 6.4
As in the Lemma 5.3 we have L(aP H.)) = VEVT
d 3 (t)
+ VEU Y = 0.
(F.66)
dt t
dt
27

---

Published in Transactions on Machine Learning Research (06/2025)
By multiplying V from left-hand side, have
d²
3 d (V Yi())) + +2V))) Y((t))
(F.67)
di2
di
where used V V = In. By replacing V Yi(t with (t),,
d
3 (d),
+ +Y)((t) = 0,
(F.68)
dia
t
dt
d' (())
3 (1) t
+
Y((1) = O.
(F.69)
dt2
dt
On
This means that for any i E [n],
de (,,(1)
3 ddYi,(()
+
+o,Yia(1) =0
(F.70)
dt2
t
dt
with (0) = -v. (1, -0
Lemma F.1. For modified Bessel funclion of the first kind I,),), have
(8)
lim ++0
1.
(F.71)
8
Proof. With Taylor expansion around 0, we have
1
t
2k+1
(1)) = Σ
(F.72)
KT + 2)
k=0
where [[k+ + is Gamma function at k + 2. For > 0,
11s))
9
(F.73)
:! 1)!
k=0
Then,
[1()) 1 1 lim
(F.74)
S
0!1!
2
Proposition F.2. ODE F.70 has the following explicit expressions:
2y1.
J (tvo), of 17>00 0.
too
(t)
(F.75)
2y
T (1 if a, 20,
w--a
where J.(-) is Bessel function of the forst kind, I()+) is the modified Bessel function of the first kind.
Proof. Firstly, consider the case of > 0. For simplicity, omit the indices I and 7((u) uY,(u//o) which satisfics
2 7 n²: (u² 1) 0.
28

---

Published in Transactions on Machine Learning Research (06/2025)
When we set Vo, we have
d2; =
YYY
11
= d
(Y)(( = -Y,()) + +YY;
du
du
Vo,
du
Zi d
Y;
2 Y
Y
u
du
du
vo /
Vo
0
0,
So, we proceed the calculation
n²244 + uZ;
+ (1²
1)
= -122 Y
u
Y,
11
}
+ u
Y
+
u
Y
1
| + (n² - -))uX,
u
Vo,
Vo,
vo,
Vo.
20² Y,
u
u
+ uY,
u
+
Y.
+ Y
u
- uY
Voi
Vo,
0
Vo,
Vol
/o,
Vo
u
+
3
Y
+ oj Y,
}
=0 0.
oi
u
Vo
Vo,
The solution function (() can be expressed with J1 (n)
(-1)
Ji(u)
u 2m+1
(2m) !!! + 2)!!
m=0
where we get around
[i( = (1 + o(1)) u
2
(F.76)
near zero for the variable Remember that Y (0) = Vi. have
(0)
2y (0)
Yis(t) tvo,
Jn (t)),), for o; > 0.
(F.77)
Secondly, we address the case of o, < 0. (u) = Y (Y,((u/Voo)) and we obtain the following modified
Bessel differential equation:
22 uZ, - (u² + 1)2: =
(F.78)
This also comes from the following calculation:
dY: 1
(F.79)
du
V σ
dZ: [YY
u
= Y
+
Y
u
(F.80)
du du
V-o;
V-o
V-o,
V-0
dt dZ, d
du² du
(Y:
u
+
- o,
u
Y;
(-m,))
u1
-
2
Y,
V-o
u
a
u
Y:
V-o
u
(F.81)
Then, we continue
122 + uZz
(u² + 1)
(F.82)
=n²
u²
2 Y,
u
u Y
u
+
Y,
+ -
Y 8
21
(F.83)
oi
0i
-
- (u² + 1)uY,
(F.84)
3
11
3u
Y
Y
41
(F.85)
oi
0
V-o
3
Y
u
3
Y
)
Y,
(ymag)
(:
(V=,)2) --o)).
(F.86)
σi
29

---

Published in Transactions on Machine Learning Research (06/2025)
Similar to the case of o; > 0, have the explicit expression of Zz(u) with I1 (2)
20
2m+1
[1(w)
Σ m=0
m! (m + 1)!
or,
(0) (t,
2y
I (tv-o;;), for 02 O.
(F.87)
W-
This formula can be extended under o = 0 with the formula F.71 i... lim [1u//uu = y
As in Lemma 6.3
0i are R-nonpositive eigenvalues, and there are (n R)-positive eigenvalues. So, Yi(1) we finally obtain
was calculated
V Yi()), and multiplying
(0) 2
2
Y,(t) = V
U1.1
I (tv 02).
I ((tv-oR)
10-02
--OOR
(0)
(0)
2+1
2y
Ji
J J(((von)
(F.88)
WOR+
W
= (u uR) 1,
2
11(tt
02
2
-OR)
(0)
vi.R)
(0)
(F.89)
1V-02
---OR
+ (uR+1...... un diag
2
Ji
2
J1 [1(t (n.)
(11, ................................ Vi.
(0)
(F.90)
VO++
Wvon
-UT (0) + -- U2T2(1Dii (1Ni,,+mm
(0)
(F.91)
where U = (uI...., E Onn, R) and U_ (uR+1...., E O(n,n n R) is orthogonal complement of U. Also
2
2
(( diag (1,
/ /1
02)
(IV--OR
-OR
(F.92)
1V-02
W-on
(2 diag
J((tVoR+1).
(1(t/on))
(F.93)
tvon
It proves the proposition.
G Implicit regularization
G.1 Proof of Proposition 6.5 concerning MM
The proof for the MM can be directly adapted from proof of Theorem 10 Cai & Ma (2022). requiring only the substitution of e -t(oA;- __1) for e 1-mm (a); - (1))
G.2 Proof of Proposition
6.5 concerning NAG
Before proceeding with the proof, let first establish the following lemma:
Lemma G.1.
1. For orthonormal basis U E O(n, R), we can define U. E O(n,n R) and oblain
|U|| = U|||=
(G.94)
2. We have the following asymptotic approrimation with J1(s) for > 1,
j)((s)
2
cos (s 3m
(G.95)
333
30

---

Published in Transactions on Machine Learning Research (06/2025)
3. Similarly, for modificd Bessel function [s)), have the asymptotic approrimation for 1,
(s)
e
(G.96)
V
Proof.
For G.94 we obtain the results with basic linear algebra.
The asymptotic expansion of (s) is found in 7.21 of Watson 1922) or 10.7.8 in Olver Maximon (2010) such as
2
3π
Ji(s)
cos (s
(G.97)
4
Dividing by O, we have the formula.
The statement follows from 7.23 of Watson 19222
The basic strategy on the proof for NAG still aligns to Theorem 10 in Cai & Ma (022) Let Uo E O(n R)
be the matrix whose columns span the null space of L(P"), and the first column of Uo is n 11/11. Also, let
E Onn, R) be the collection of eigenvectors of L(P) corresponding the smallest R eigenvalues. With the Davis-Kahan theorem discussed in von Luxburg (007
Yu al. (2015), we have
[L(E]]
1 U|l = i A
(G.98)
where E= -P-Pa and Uoi is orthogonal complement of Up.
U.UJ (1) Y(OIIla
(G.99)
(Y)(0)
10.U0 UT (t)yt)). (0) UT UT (tly O0)
UJT ( U.F2(t)y(RR (0)
ll
(G.100)
|Y((0)]
11U..U - (OUTT(Uy)L..))
[U.U I)U
(2(0y L..ell
(G.101)
||Yx())lle
(0Y
11(101 U._ JUT (
11 U. J
(2
(0)
(G.102)
(Y)()l
||Yl())
||U
(t)l
uT
(G.103)
| I
(G.104)
||(EE]| 2
2
11 (1(V-ORR) +
J (1(W/R+1)
(G.105)
(L+)(L(P))) W-OR
WOR+I
|L(E]]
2
11 (aXR(((P))) +
11
1
(G.106)
(LPx)((P)) /1 -add((((P)) +
2 J (1 (L(P)) 1
(G.107)
W (LPP) -
-
Whenever
[[(E)]] & AL+ (a)n(((())) 1 --1
(G.108)
31

---

Published in Transactions on Machine Learning Research (06/2025)
with ())+)(L(P)) > (fdx+11(LLPP) + Jo, we have
U Yi(() Y1(O)22 0.
lim (t,n)-oc
(G.109)
(1(0)]2z
Also,
[UGU Y(()) - U (,)l2 Y(0)l
(G.110)
Y2(0)]l2 10.00 UT (t)y U. (0)
U.U, UT U - U,
la
Y/(0)]22
(G.111)
10.00 (() U(T (())))) (0)
10.0 U (T) (2 (0)
(G.112)
|Y)(0)12
Y((0))22
IT(() (t) (0))]l + I U (T2(t) (2(0))
(G.113)
2
[[(E]]
2
1i 11(1V-an
Ji (I(1/OR+))
(G.114)
W on
AR+1 (L(())
Wont
2
11 (1(4 -adn(L(P))
n n-1 n-1)
1
(G.115)
ann (LP())
n
|L(E]|
2
J (()))+I(L(PP)
(G.116)
AR+1(((P')
t I
n-1
Whenever the condition G.108) holds, we also have
-, Y,00 0.
lim (t.n)
(G.117)
|Y,(0)]12
Eq. G.108) is ensured by the conditions of the theorem and the asymptotic behavior in F.71 G.95
follows:
The condition [[(E)]
(fx [L(P' P)[I} & (cA)(L(P))- + comes from LPP -P]] as n + x in (T2.N) and t t= o(n 12) in (TINN).
(L+1(L(P))) comes from AR+1 (L((")))))
Also, with the Weyl's incquality I):(L(P)) (L(())) < [[EE]]| (I 2 aL(P") | > 1
due to >> (nAa+1(L(P")) in (TI.N), Finally, t2 f + co can be deduced with Weyl's inequality and A (L(()) >> (Po) in (T2.N).
Lyapunov exponents for each eigensolution
The expressions
16 or 24 can be interpreted as the temporal evolution of a dynamical system. By
examining the Lyapunov exponents, one can gain insight into the asymptotic bohaviors for
H.1 Lyapunov exponent for GD and MM
We calculate the Lyapunov exponent for MM We pick up the -th-term from the equation 16 and set Lyapunov exponent as Atyapa for E [n]:
J (2,((1
ALyap,
lim
lim
log
(H.118)
1-+00
8Y,
(0)12
For GD, consistently interpret m throughout this subsection.
32

---

Published in Transactions on Machine Learning Research (06/2025)
where (t) the perturbation at time t caused by the initial perturbation oY,,, Now, we have
(i,(t) exp
ad;
n2 )
1
(0) uu yi Ju, [2].
(H.119)
m
Therefore, by canceling out the same terms in the numerator and the denominator,
lim lim Aayap, (-++0 (0) +0 log exp ( (-1
1 m n -111)
1
aa:
(H.120)
1 (a1;
(H.121)
Note that the result holds regardless of E [2). With the Lemma 5.5 have the following proposition Proposition H.1. (Iyaauuoo exponent. for GD and MM) Under the conditions (T1.M), (T2.M), >> and i E [n].
m (ax Especially, ALyop. 20 (i 1,--,,.. 1,..., ... R), and Aggon. < 0 (i = R + I....,n. n).
Mayan'
(H.122)
This statement suggests that the elements of Yi(() (i l,..., R) highlight a specific structure within the cluster, enhancing its distinctive characteristics.
H.2 Lyapunov exponent for NAG
Now, for the case of NAG. Using the formula (24 the asymptotic approximation G.95 and G.96). we arrive at the following conclusion:
Proposition H.2. (Lyapunov exponent for NAG) Under the conditions (T1NN), (T2 and > 1, have
add+ + n n- (i= R)....
ALyapá
(H.123)
0,
(i R m))
Especially, ALyap, 2 Of1 = 1,..., R).
It also suggests that the elements of (t) (i (i= .,....,I, Rt) emphasize a cluster structure.
The case for R+I....,, .n differs slightly. GD or MM, ALypp, is negative, causing to exponential decay. However, in NAG, oscillations occurs, the Lyapunov exponent is considered to be 0. This aligns with the typical differences between GD, MM and NAG.
Supplemental numerical experiments
1.1 Differences between solution paths of ODE corresponding to GD, MM and NAG
We dorived ODEs corrosponding to the iterative optimization procedures for t-SNE with and without accel- eration. We also derived their closed-form solution. We compare the difference between solution paths by simply plotting them in Fig. 4 As illustrated Fig for E [R], i.c., components corresponding to the negative eigenvalues of the Laplacian L(aP H.) (let's denote o < 0), GD follows an exponential function exp(-to), and the MM follows exp(-1=t et
- both of which are exponential functions (see Eq. (16)). On the
other hand, NAG follows function [/ -
4/-o involving the modified Bessel function of the first kind (1(-)
(see Eq. 24 Based on their series expansions, these functions are greater than for positive values of and converge to as approaches zero (see Lemma for NAG case).
By contrast, Fig. for > 1, GD and MM docay exponentially and NAG decays under the for- mula G.95). From these conditions, the cigenvectors u, for 2 R 1, ie., components corresponding positive eigenvalues, asymptotically decay to as time increases.
33

---

Published in Transactions on Machine Learning Research (06/2025)
1.0
0.8
30
9.6
Fuanteen
1.003
Jundaan
0.2
0.0-
1000
(a) GD, MM and NAG (negative eigenvalues)
(b) GD, MM and NAG (positive eigenvalues)
Figure 4: Comparison of (Modified) Bessel function and exponential functions.
1.2 Experiments with synthetic dataset
We conduct a comparison between iterative optimization methods for t-SNE and ODEs corresponding to them, specifically for GD, MM, and NAG, to explore their convergence behaviors using a synthetic dataset, generated from the following Gaussian Mixture Model (GMM):
X T (u,, Σ), Zi Multinomial T2, T3), for i [200]:
(I.124)
where μ1 (0, f2
(150, 110, 170) 43 = (-130, 1 -150) are the centers of Gaussian distri-
330
20
25
butions, =
20
50
10
is the common covariance matrix and (TT1, m2, R3) =
With the
25
10
30/
above (200)_ we define the element of P as Pij = Ply if 2 = 2y: otherwise, is set 0. The Laplacian
L(oP H has an eigenvalue with almost and other two negative eigenvalues, i.c. R 3. According
to Cai & Ma (2022), it can be proven that the conditions (TI.M), (T2.M), (TIN) and (T2.N) are satisfied
by placing additional assumptions on these.
Figure shows that the comparison between iterative optimization methods with and without acceleration (GD, MM, and NAG), and ODEs correspond to them for the dataset generated from the above explained GMM. The optimization done with step size parameter
5, momentum parameter m 0.5, Perplexity is 30, and exaggeration parameter a 10. In the plot, solid lines indicate iterative approaches and dashed lines denote methods using ODEs. From this plot, the I-SNE optimization procedure) aligns well wetw iterative algorithms and continuous limit ODEs. It can also be observed that as increases, the discrepancy in the KL-divergence values between those obtained from the iterative algorithm and those obtained as the solution path of the ODE becomes larger. This is consistent with the theoretical results obtained in this paper.
that the KI--divergence (the objective function of
As in section we state the assumption Lemma 5.1 about GD and MM, we identify the variables as kh, where the variable & as time parameter, k iteration number and h step size. For NAG, the identification can be t kvh as in Lemma 6.1 Figure 6 shows the representation where the top row corresponds to GD, the middle row to MM, and the bottom to NAG.
We have two kinds of observations.
Case 1: Focus on that the middle case in GD (t = 35) and the rightmost case in NAG (t = 33.54). In the case of NAG, even though the value of t is smaller at 33.54 compared to 35 in the case of GD, the data well-separated, whereas in GD is still undifferentiated..
Case 2:
The rightmost case in GD (L = 75) and the middle case in MM (1 (1= 35). Formally, the time parameter is different, and the dineeninsio representation looks similar. Considering the effect of the momentum term m 0.5, can be interpreted as accelerating time t, effectively making 70(= 35/(1- (1)5)) See the formula (66), where we adopt that MM can be understood fast-forwarding GD with respect to the variable thanks to the moment term with m.
34

---

Published in Transactions on Machine Learning Research (06/2025)
12
14-
KL-divergencs gevergencc
0.8
0.6
GD_ODE ODD
GD iterative
MM_ODE
iterative
NAG
NAG
60
80
100
Figure 5: A comparison between iterative optimization methods and those approximated using ODEs (GD,
MM, NAG).
1.3 Profiles of datasets
addition of GMM datasets explained in section 1.2 the following datasets are used. Here are the links
and their licenses:
KDDCup1999: cense referenced
http::////aaaaaaa..////////////a CC BY 4.0 li-
MNIST: http:///stttteoeee..../..///a//s
openm1 .html BSD license provided by sklearn
Olivelti Faces:
fetch_ brrccholiveetiifeeee faces. html
BSD license provided by sklearn
For the KDDCup1999 1999 dataset, download the keddup. dala 10 percent.gz file locally and it. For the MNIST and Olivetti Faces datasets, download and use them directly in the code by utilizing the functions provided by sklearn.
1.4 Experiments through variation of random Initialization
It is well-known that the clustering results using t-SNE may vary with variations of initialization such as Kobak & Linderman (021)) Cai & Ma 2022) states sufficient condition in their Theorem 14 as random initialization which leads to intercluster repulsion (Th.13) after both early exaggeration and embedding stage. It is reported that false clustering may appear due to an incidental combination of overlapped cluster from the early exaggeration stage as in Remark 16 there. By contrast, we focus only on early exaggeration stage, and numerical experiments on several datasets, performed multiple runs with different initializations for the same dataset, and present the evaluation results.
We use KDDCup1999, MNIST and Olivetti datasets. The Olivelti face data consists of images of 40 individuals with small variations in viewpoint. The data set consists of 400 images (10 per individual) of size 64 x 64 = 4096 and is labeled according to identity.
In order to evaluate the variation of initialization, we used Adjusted Rand Index (ARI) (Hubert & Arabie 1985), which commonly employed to evaluate clustering performance with known labels (ground truth clusters). This index assesses the degree of correspondence between the predicted clusters and the true class labels, serving as a measure of clustering accuracy cases where the true labels are available, akin supervised learning.
35

---

Published in Transactions on Machine Learning Research (06/2025)
GD (k=4.h=54 5,0
GD (km 5,
GD (k= 15,-15,hh5555 75)
0.06
0.04
0.02
0.00
0.02
-0.04
-0.050
-0.025
0.000
0.025
-0.02
0,00
0.02
-0.01
0.00
0.01
MM hk=4.h=5,t=20)
MM (k=7.h=5.t=35)
MM (k=15
0.04
0.03
0.02
0.01
0.00
-0.01
-0.02
-0.02
0.00
0.02
00.11
0.00
0.01
00.1
0.00
0.01
NAG (k=
NAG (5566)))
NAG (k=155=55
0.100
0.075
0.050
0.025
0.000
-00055
-0.050
-0.075
-0.10
-0.05
0.00
0.05
00.04
-0.02
0.00
0.02
-0.02
0.00
0.02
Figure 6: Low-dimensional expression with GMM dataset. Top: GD, Middle: MM, Boltom: NAG.
Figure 7 demonstrates the results. The band visible in the ARI transition represents the results of varying ARR after randomly changing the initial values 30 times. The width of the band corresponds to the standard deviation of the ARR.
In all three datasets, ARI increases most rapidly for NAG(green), followed by MM(purple), with GD(pink) trailing behind. NAG achieves the best results the shortest time. It important to note the relationship between the variables and for GD and MM, kh, while for NAG, = kvh. As result, even with the same value of k the effective value of & is smaller for NAG, which makes the NAG curve appear truncated.
36

---

Published in Transactions on Machine Learning Research (06/2025)
(a) KDDCup1999
(b) MNIST
(c) Olivetti
0.8
0.6
GD
033
0.3-
0.6-
0.4
Atture Kans
0.4
03
0.2
GD
400
600
800
200
400
3K
340
50c
Figure 7: Adjusted Rand (Ade for GD, MM, NAG across different dalasets
1.5 Experiments through various initialization methods beyond random initialization
Generally, there are several ways to initializing methods for I-SNE. Kobak & Linderman (2021) reports the result with initialization using PCA, Linderman & Steinerberger (2019) uses spectral clustering in initializing. In this section, we experiment several initialization methods and compare the results under the equation 16 or (24] with the indicator ARI.
We basically follow the way of section 1.4 The MNIST dataset is used for the experiments, and the Adjusted Rand Index (ART) is employod as the evaluation metric. To assess the stability of each initialization method, conduct 10 trials for each approach. In addition to random initialization as described in section I.4 we attempt. representative initializations using PCA, Spectral Embedding (SE) (Belkin & Niyogi 2003). and
MultiiDmennionall Scaling (MDS) [Torgerson] 1952p. For PCA, SE, and MDS, the dimension is reduced to two by setting n components to 2. In the case of SE, the number of neighbors (n (a_neighbors) used to construct the Laplacian matrix is set 300.
Figure 3 demonstrate the results. For the optimization methods GD(orange) and MM(purple), we observe that changing the initialization from random other methods results in starting with relatively higher ARI, and the metric generally improved the time increase.
On the other hand, for NAG (green), it is observed that the ARI starts near zero in all initialization methods. This is because, in the case of GD and MM, all components of the initial values are included in the eigendecomposstionn of Eq. 16) like u,Tr1),
(0) whereas for NAG, the decomposition in Eq. (24 only includes In the case of NAG, as time progresses, the influence of the eigenvectors becomes reflected, and the ARI gradually improves.
a subset of the initial values in each term like ( or U_T2(t)y(
Furthermore, across the different initialization methods, is confirmed that the SE leads the greatest improvement in terms of ART in experiment.
As minor point, for MDS, using random placement during initialization results in similar variation in ARI outcomes as observed with random initialization. In contrast, PCA and SE, being deterministic methods, do not exhibit such variation.
1.6 Experiments through ARR for other data sets
We also explore the stopping time obtain the clustering results with ARR(Average Residual Rate) defined in (88 at appropriate time. In addition, we observe whether another indicator 'Trustworthiness exhibits a similar trend to ARR during the optimization process.
In section
evaluate ARR for MNIST data set. Here, we evaluate it with KDDCup1999, Olivetti Face datasets as well. Figure g shows the results for KDDCup1999, and Fig. 10 for Olivetti.
8.2
7 "Remember that both (1) and (2( are diagonal matrices Theorem 6.4
37

---

Published in Transactions on Machine Learning Research (06/2025)
random
optimization methed
GD
36-
NAG
Adjusted
9--
43-
optimication method
A2+
MM
NAG
293
413
02.
80C
50
mds
Adjusted Rand (ARI)
0.4
sptimication method
optimication method
GD
GD
NAG
A1-
NAG
42
Figure 8: A comparison among various initialization methods with MNIST dataset for GD/MM/NAG
The left figure shows the transition of ARR and Trustworthiness. We begin by reviewing ARR. (AG(green) exhibits the fastest drop in ARR, followed by MM(orange) and then GD(blue). A red dashed line marks the threshold at ARR 0.01. As mentioned in section L.4 it is important to note that due to the different relationships between the variables and & among GD M the NAG curve appears truncated.
Next, we compare ARR with Trustworthiness (Venna & Kaski 2001). Trustworthiness is a standard metric used in methods such as t-SNE to evaluate how well the structure of the original high-dimensional data is preserved in the low-dimensional embedding. Specifically, it measures whether the neighbors the embedded space are also neighbors in the original high-dimensional space. The value ranges from to 1, with values closer to indicating higher agreement.
The same figure also shows the case where the number of neighbors is set to 10 for Trustworthiness. The dashed GD/MM/NAG. As ARR decreases, Trustworthiness tends to increase, indicating an approximate inverse correlation. The correlation coefficient around -0.85 for KDDCup1999 and -0,90 for Olivetti. Further- more, because the Trustworthiness value is sensitive to the choice of neighborhood size, the strong inverse correlation observed in our results does not necessarily indicate that trustworthiness can be used as general alternative to ARR. for determining the stopping time.
green lines represent the Trustworthiness values for each optimization method:
Although Continuity is also a well-known metric (Kski et al. 2003) similar to Trustworthiness, numerical experiments in this case showed that its values were nearly equovalenttthough not identical to those of Trustworthiness. Therefore, we only present Trustworthiness in the figure.
The right panels of figures and 10 show the clustering results of GD, MM, and NAG when ARR - 0.01.
1.7 Experiments with various momentum coefficients
We explore the effect of the momentum coefficient m in the MM. In the Eq. 16 the transition of ARR for MNIST dataset when changing the coefficient m from 0.1 to 0.9 can be seen in Fig. Other configurations of MNIST are the same as in section 8.2 The cocfficient m was originally assumed to range from 0 to 1, and can be confirmed that as m increases, the influence of the residual term on ARR diminishes over a shorter period of time. Formally, GD can be regarded as having m of 0, which is also consistent with the results shown with both datasets -KDDCup199 and Olivetti- in Figs. g 10
8 kh for GD and MM, whereas kvT for NAG.
38

---

Published in Transactions on Machine Learning Research (06/2025)
ARR Wisshwortsissssstt for KDDCup1903 dirsset
18
(1=500)
MWit-250
(t-10436)
0.0
ARN Trusttorthhi wwrssttttiinsss
Method Netii
CD TW
MM TW
NAG TW
112
-00 600
-004 460
Figure 9: (left): transition ARR and Trustworthiness (TW), (right): Clustering result for KDDCup1999 dataset with ARR=0001 using three ODEs: GD, MM, and NAG.
(b=1030)
MMLit=520)
NAS (t-192.9)
ANR TTugwwrrthness
Method Metric
CD ARR
6.0
GD
MM ARR
MM
NAG ARR
2008 425
020 COCERIS (20) SOW 3000
8.8 cioc aa1o 2015 1ies R021
Figure 10: (left): transition of ARR and Trustworthiness (TW), (right): Clustering result for Olivetti Face dataset with ARR=0.01 using three ODEs: GD, MM, and NAG.
1.8 Further analysis of KDDCup1999 dataset
We obtain the adjacency matrix P as depicted in Fig. 12 The heatmap illustrates that the elements of P exhibit relatively high values (orange or yellow) between clusters, while other regions are closer to black.
Solving the eigenproblem with the Laplacian matrix L(aP Hn) using a = 10, we obtain eigenvalues close zero and four negative cigenvalues, as shown in Fig. They are emphasized during clustering regardless of the method used (GD/MM/NAG), resulting in the formation of clusters in the low-dimensional representation.
Furthermore, Fig 14 displays the distribution of the dominant eigenvectors U1,..., These cigenvectors are emphasized in the low-dimensional representation as the time parameter progresses under Eqs. and 24 As shown in this figure, the components belonging to the same cluster simultaneously take on the same values, and the relative magnitudes of these values become crucial when forming clusters in the low- dimensional space. This mechanism leads the clustering results depicted Fig. the oiv-dimensional representation.
39

---

Published in Transactions on Machine Learning Research (06/2025)
with MM MNIST dataset
0.5
momentumn coefficient
0.1
Avarage Ratio(AR))
0.3
0.4
0.7
0.9
0.2
0
500
1000
1500
2000
2500
3000
Figure 11: A comparison of MM with various momentum coefficients
100
300
2
-1
100
200
300
400
Figure 12: A heatmap of similarity matrix P for the 500 samples from KDD Cup 1999 dataset.
40

---

Published in Transactions on Machine Learning Research (06/2025)
values
values of
values of
0.025
Distribution of Eigenvalues
0.000
-0 -0.00
Count
0.075
40
0.050
0.22
0..50
-0.75
-0,002
0.000
0.002
0.004
0.006
0.008
0.010
0.012
200
Figure Distribution of eigenvalues (aiheen)) of the Laplacian matrix L(aP H,): of O, of 0(2 i 5) and > 0(i 6).
Figure 14: The eigenvectors, whose eigenvalues are the smallest five: u1,..., us.
41


Reviews:
[{"review_id": "m9kqna7OGb", "reviewer": "Reviewer_qjTq", "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\n1. the theoretical derivations, motivation and contributions are in depth and non-trivial, although it can benefit from better exposition. \n2. The empirical work gives new insights and is relevant to the corresponding theory\n3. The topic is relevant, important and potentially of broad interest \n\nWeaknesses. \n\n1. The exposition is dense and lack enough motivation and big picture.\n2. Unclear what part of the analysis is novel versus modifications of prior work.", "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": "The paper investigates the dynamics of t-SNE. It advances our understanding of t-SNE using the perspective of viewing iterative optimization as ODEs and investigate the convergence behavior under accelerated gradients/momentum.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- Consider giving some type of outline or big picture diagram in the appendix or intro that gives a summary of how the components fit together (i.e. what you are trying to show, what tools allow you to show that, and why those tools are appropriate etc) \n\n- Consider adding a small paragraph outlining exactly which technical components are new and what are borrowed from prior work so it is clear to the reader \n\nMinor changes: \n\n1. Since this is an ML audience, use of physics terminology such as ansatz might not be universally understood in the community. Perhaps can add some clarification/definition here to make it easier for the ML reader, and to clarify what exactly the notation $\\approx$ here is referring to.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "LwolgcUJhc", "reviewer": "Reviewer_xx4b", "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  \n1. The paper develops nontrivial connections between gradient-based optimization and continuous dynamical systems involving Bessel and modified Bessel functions. The derivation of the ODEs under MM and NAG is rigorous.  \n2. The approach allows analytical computation of embeddings at any time 𝑡, which is conceptually elegant.  \n3.  The proposed Average Residual Ratio (ARR) is a clever and interpretable metric tied to spectral concentration.  \n\nWeaknesses  \n1. The analysis closely parallels to [1], only extending from GD to MM/NAG. While nontrivial technically, is this sufficient for a full TMLR paper? The key observation that acceleration introduces Bessel-type behavior is interesting but lacks broader implications or surprising insight into t-SNE dynamics. \n2. The use of eigendecomposition of 𝐿(𝑃) implies $O(N^3)$ complexity, which contradicts the paper's emphasis on acceleration and efficiency. How does this framework apply to realistic high-dimensional datasets like CIFAR or single-cell RNA-seq with 10K–1M samples?\n3. All experiments are limited to small datasets (MNIST, KDDCup1999, 1600 samples). These are too simple for demonstrating any practical edge of the proposed method. \n4. The entire analysis is confined to the early exaggeration (EE) stage. However, meaningful embeddings often require a post-EE embedding phase. The paper would benefit from analyzing or adapting the embedding phase itself.  \n\n[1] T. Tony Cai and Rong Ma. Theoretical foundations of t-SNE for visualizing high-dimensional clustered data.", "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 presents a theoretically motivated extension of the dynamical systems view of the t-SNE algorithm by incorporating momentum (MM) and Nesterov’s Accelerated Gradient (NAG). The authors derive closed-form ODE solutions, reveal connections to Bessel functions, and propose an early stopping criterion (ARR) based on spectral decomposition. The approach is novel in framing accelerated t-SNE through analytically solvable continuous-time dynamics, and the paper is well-supported by derivations and empirical results.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "1. What new understanding of t-SNE's clustering behavior does the proposed method bring that is not already accessible via the GD-based analysis?  \n\n2. How does the proposed method handle complexity problem?\n\n3. Can you demonstrate a real-world case where MM or NAG-based t-SNE materially outperforms both standard t-SNE and UMAP?", "broader_impact_concerns": "I do not foresee any potential major ethical problem.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "1IZICem2Xj", "reviewer": "Reviewer_8vBG", "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:\n1) The authors provide a theoretical analysis of a widely used method.\n2) The closed-form ODE solutions should be of practical benefit.\n\nWeaknesses:\n1) I found the paper very difficult to follow. This is for one certainly because I am not familiar with t-SNE. But it is also due to the writing, which is very dense, does not provide much intuition, and is missing some definitions. Also the language should be improved.\n2) The paper is missing a problem formulation.\n3) It is not always exactly clear what is taken from related literature and what is a novel contribution.", "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": "The paper analyzes the t-stochastic neighbor embedding (t-SNE) mwrhos. The key contribution is the analysis of t-SNE with included momentum term from the dynamical systems perspective. The analysis shows better convergence rates of the variant with momentum term. The authors further extend to applying Nesterov's accelerated gradient (NAG) method to t-SNE to improve convergence speed further.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "1) In the literature review section, it would be good to spell out how you advance the state-of-the-art. Particularly, in Sec. 2.4, how do you solve the problem you identify?\n2) The writing is very dense and focuses on the technical concepts. It would be great to provide some intuition and explain what the different concepts mean. What is the impact of assumptions. What is the impact of the results you derive.\n3) Please make sure all elements are properly introduced and defined. For instance, in the beginning of Sec. 4, I don't think H_n was introduced. Also, why is the reformulation for which you use H_n important?\n4) Can you add a concrete problem formulation?\n5) Can you clarify whether Theorem 4.1 is taken from the reference (Cai & Ma, 2022) or has some novel elements?\n6) The concept of a stopping time has a specific definition in stochastic analysis. But to my understanding, you are using it more as a general criterion for when to stop an algorithm. I think it would be good to make that clear.\n7) The language should be improved. There are some plural 's' missing, e.g., in the first and second sentence of the second paragraph in the introduction, some articles are missing, e.g., when defining the adjacency matrix after (4), \"where put m = 0 for GD\" sounds a bit off.", "broader_impact_concerns": "No concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "znuHpf7iNE", "reviewer": "Reviewer_cTTp", "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  + t-SNE is a highly, highly used procedure, and its full understanding as eluded the community (for the most part) since its introduction, so work on this topic is important\n  + The analysis is mathematically rigorous, involves blending of understanding from several domains, and ties it together in a careful way\n  + The main mathematical insights are supported and enhanced by simple experiments \n  + The use of momentum and NAG are demonstrated to improve the efficiency, and explained why and when.  \n\n\nWeaknesses: \n  + The results are explained in a highly technical way that will mostly only be accessible to experts in both t-SNE and optimization.  While this is an appropriate framing for others who would investigate these issues, it will have less impact on the much, much broader community that simply use t-SNE (mostly as a black box).  \n  + Much of the work is an extension of the work of Cai & Ma (JMLR 2022).", "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": "t-SNE is probably the most popular non-linear dimensionality reduction method within data science for visualizing high-dimensional data in 2 dimensions.  Despite this popularity, it is not understood what formally can be concluded about its resulting embeddings, and how to understand them.  This paper improves this understanding by studying its iterative dynamics -- in particular, how it is affected by using momentum or Nesterov's accelerated gradient.  \n\nThe main results show that (and how, when) these optimization improvements help -- in general they lead to faster convergence.  It also exposes (through a closed for solution of an ODE of the system) a time parameter, and how the resulting embeddings strongly depends on it.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Suggestions:  \n I think the potential impact of the paper could be *significantly* improved if the authors add a subsection (after \"main contributions\" before Sec 2), about something like \"Take-a-ways for Users of t-SNE\".  I am imagine something like:\n  - if t-SNE is too slow for your data set (and some easy to understand condition holds) run with MM  (or NAG)\n  - The time parameter (explain in very simple terms, and how one can tune it in software) plays an important role, similar to the bandwidth choice.  Set it large for tighter clusters [or some similar broad guidance].  \n\nI trust the authors to improve upon my sketched ideas above.", "broader_impact_concerns": "none", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The manuscript has received generally positive evaluations, highlighting significant theoretical contributions and clear connections established between accelerated gradient methods and dynamical systems theory. The reviewers have confirmed that most of their concerns have been adequately addressed. 

However, as pointed out by one reviewer, there remains an important aspect that requires additional clarification.
Comparison of ARR with standard quality metrics: The authors have partially addressed the requested comparison between the proposed Average Residual Ratio (ARR) metric and established quality metrics commonly used in evaluating t-SNE embeddings. Although the current version includes numerical experiments and visualization (Figures 9, 10, 11), it lacks a direct and explicit comparative discussion or quantitative analysis with standard metrics. For the final submission, the authors must explicitly address this concern by either providing a direct comparison or clearly articulating how ARR relates quantitatively and qualitatively to existing standard metrics. This will substantially enhance the practical relevance and interpretability of the ARR metric for practitioners.
Decision: Accept (poster)