Scalable Second-order Riemannian Optimization for $K$-means Clustering
Authors: Peng Xu
Chun Ying Hou
Xiaohui Chen
Richard Y. Zhang
Venue: ICLR 2026 Poster (2026)
Keywords: K-means clustering
manifold optimization
Newton's method
nonconvex
Primary area: optimization
TL;DR: Smooth Riemannian formulation for $K$-means clustering; rapid convergence to second-order critical points with linear per-iteration cost.

Abstract: Clustering is a hard discrete optimization problem. Nonconvex approaches such as low-rank semidefinite programming (SDP) have recently demonstrated promising statistical and local algorithmic guarantees for cluster recovery. Due to the combinatorial structure of the $K$-means clustering problem, current relaxation algorithms struggle to balance their constraint feasibility and objective optimality, presenting tremendous challenges in computing the second-order critical points with rigorous guarantees. In this paper, we provide a new formulation of the $K$-means problem as a smooth unconstrained optimization over a submanifold and characterize its Riemannian structures to allow it to be solved using a second-order cubic-regularized Riemannian Newton algorithm. By factorizing the $K$-means manifold into a product manifold, we show how each Newton subproblem can be solved in linear time. Our numerical experiments show that the proposed method converges significantly faster than the state-of-the-art first-order nonnegative low-rank factorization method, while achieving similarly optimal statistical accuracy.

Published as a conference paper at ICLR 2026
SCALABLE SECOND-ORDER RIEMANNIAN OPTIMIZA- TION FOR K-MAANS CLUSTERING
Peng Xu
Chun-Ying Hou"
University of Illinois Urbana-Champaiga pengrul@illinois
University of Illinois Urbana-Champaign
mroouu@illinoi....
Xiaohui Chen
Richard Y. ZZhang
University of Southern California
University of Illinois Urbana- Shrana-Cammpaig
xiaohuic@use
ryz@illinois.ee
ABSTRACT
Clustering is hard discrete optimization problem. Nonconvex approaches such low-rank semidefinite programming (SDP) have recently demonstrated promising statistical and local algorithmic guarantees for cluster recovery. Due the combi- natorial structure of the K K-means clustering problem, current relaxation algorithms struggle balance their constraint feasibility and objective optimality, presenting tremendous challenges in computing the second-order critical points with rigorous guarantees. In this paper, we provide new formulation of the K K-means prob- lem as a smooth unconstrained optimization over a submanifold and characterize Riemannian structures allow to solved using a second-order cubic- regularized Riemannian Newton algorithm. By factorizing the K-means manifold into a product manifold, we show how each Newton subproblem can be solved in linear time. Our numerical experiments show that the proposed method con- verges significantly faster than the state-of-the-ar first-order nonnegative low-rank factorization method, while achieving similarly optimal statistical accuracy.
INTRODUCTION
Clustering is cornerstone of modern unsupervised learning, where the goal is to group similar observations into meaningful clusters. The problem commonly approached through the K-means formulation, which seeks to partition n data points X1, X2...., X, E IRd into K disjoint groups G GI...., GK by maximizing the total intra-cluster similarity:
max G.....GKK
Σ IJEGG
(X X3)
G [n]
GR
Here, the inner product X,) X,Tx, is used to measure pairwise similarity, G denotes the cardinality of Gk. and Li denotes disjoint union. Most common algorithms for K K-means clustering, including Lloyd's algorithm (Lloyd, 1982) and spectral clustering (Ng et al., 2001 von Luxburg 2007). can be understood as heuristics for finding "good enough" solutions to the discrete optimization
(1) These methods do not come with any guarantees of local optimality, let alone global optimality. Indeed, it is commonly argued that globally solving (1) is NP-hard in the worst-case (Dasgupta) 2007
Aloise et al. 2009), and would lead to statistically meaningless clustering that overfts the data.
Yet in average-case regimes, globally solving the K-means optimization problem computationally tractable as well as statistically optimal. In particular, when the data X,,..., X,, arise from Gaussian mixture model with sufficiently well-separated components, Chen Yang (202bb) showed that well-known semidefinite programming (SDP) relaxation of Peng & (2007
be both
"These authors contributed equally.

---

Published as a conference paper at ICLR 2026
written
(XXT, 22 + Σ max ZERuXn
log( 2), 21, = =1, tr(Z)= K, Z -=
where == (X1,... XT) and Zz.) + = log(max (Zij: 0)), ), guaranteed to compute the globally optimal clusters Gi...., G for in the limit - 0+ in polynomial time, that in turn recover the ground truth partitions. Note that the formulation (2) is equivalent to the standard K- means SDP formulation with the elementwise nonnegativity constraint Z 2 0 in Peng & Wei (2007) Chen & Yang 2021b) (see Appendix A for more discussions). Moreover, this recovery occurs
as the separation between the clusters is large enough for it be possible. Put in another way, if solving (2) does not recover the ground truth partitions, then the clusters are too closely spaced in way that makes recovery inherently impossible in information--theortticc limit sense, see Section 2.1 for more details.
Unfortunately, the SDP ( is not practical means of solving (1) to global optimality, due its need to optimize over an X matrix to cluster n samples. Following Burer & Monteiro (2003); Boumal al. (2020). natural alternative is factor UUT into its r factor matrix U for rank parameter r 2 K, impose the logarithmic penalty over U instead of Z, and then directly optimize
over U:
(XX ,UUT) + + [og(U : UU In tr(UUT) K }
max UEROXP
This reduces the number of variables and constraints from O(n²) down to (n)). but at the cost of giving up the convexity of the SDP. In general, we can best hope compute critical points, which may be spurious local minima or saddle points. The core motivation for our approach, and the impetus for this paper, is the surprising empirical observation that all second-order critical points are global optima in this setting; this is formalized as the following assumption.
Assumption (Benign nonconvexity). In the average-case regime when (2) globally solves ( all approximate second-order critical points in (3) are within neighborhood of a global optimum.
The phenomenon of benign nonconvexity is well-documented in the unconstrained version- optimizing over semidefinite Z 0 by factorizing Z UU datingg back the early works of Burer& & Monteiro (2003). In contrast, it is rarely seen in our nomegative variant, which adds the elementwise constraint U 0 enforce doubly nonnegativity in Z UUT, Despite superficial similarity, the two formulations differ in fundamental ways, with the nonnegative case known to admit numerous spurious critical points; see Section 1.2 for some classic and recent examples. Nevertheless, we consistently observe that all second-order critical points correspond to global optima, that in turn successfully recover the optimal clusters.
1.1 CONTRIBUTIONS: CHEAP AND FAST CONVERGENCE TO SECOND-ORDER CRITICAL POINTS
Under Assumption 1 globally solving the K-means optimization problem ( reduces to that of computing second- critical point for (3). Unfortunately, the constrained nonconvex setting, there is no general-purpost algorithm that is rigorously guaranteed to compute a second-order critical point. The core issue the need maintain feasibility, i.e. for each iterate U to satisfy the K, while making progress towards optimality.
nonconvex constraints UU 1, = 1, and trU General purposs solvers like fmincon (Byrd et al. 2000) and kni (Byrd al.) (2066) promise convergence only to critical points of an underlying merit function, which may be infeasible for the original problem. Augmented I Lagrangian methods guarantee convergence only to first-order critical points, and only when starting within a local neighborhood (Zhuang et al. 2024). This is significant departure from the unconstrained nonconvex setting, where a diverse range of algorithms both cheap first-order algorithms like gradient descent, well aapilly-coverre trust-region Newton's 's method globally converge to a second-orde -odder critical point starting from any initial point.
second-order methods like
Our first contribution is to present an interpretation of (3) as a smooth unconstrained optimization Riemannian manifold. This allows the immediate benefit of extending the wide array of
2

---

Published as a conference paper at ICLR 2026
unconstrained optimization algorithms to the constrained setting, well as their accompanying guarantees for first- and second-orded optimality. For the first time in the context of K K-means, we open the possibilitty to guarantee global convergence to first- and second-order optimality.
Our second contribution is to show that second-order Riemannian algorithms can be implemented with linear per- per-iteration costs with respect to the number of samples n. In other words, of all practical algorithms to compute second-order critical points, we show that the one with the best iteration complexity (second-order methods) can be improved to have the same per-iteration costs as first-order methods. Our final algorithm computes second-order points in e -3/22, poly(r, d) time.
1.2 RELATED WORK
Benign nonconvexity in the unconstrained Burer-Monteiro factorization Z - UUT has been empirically observed since the early 2000s (Burer & Monteiro 2003), and widely exploited in nonconvex low-rank algorithms in machine learning. In the past decade, theory has been developed to
explain this phenomenon under some specialized settings (Bhojanapalli et al.) 2016 Ge et al.) 2016 Bandeira et al. 2016; Boumal et al.) 2016 Ge et al. (017). Unfortunately, these guarantees tend to be conservative in the number of samples or the level of noise; they capture the general phenomenon but cannot rigorously explain what is broadly observed in practice.
In contrast, the nonnegative Burer-Monteiro factorization = UUT with U 2 0 is widely understood not exhibit benign nonconvexity. To give two simple examples, the functions are easily confirmed to exhibit benign non- convexity over U E Rnxrr But imposing 2 0 causes spurious local minima proliferate; this is unsurprising because both problems, namely copositive testing (Murty & Kabadi 1987) and complete positive testing (Dickinson Gijben 2014). are well-known to be NP-hard. For a more waeeessoopstittte sophisticated example, the function ((0) AATUU U,U ||2 is well known to exhibit benign nonconvexity when the linear operator A: : 55t(((: + R satisfies the restricted isometry property (RIP) (Bhojana- palli et al. 2016). In this context, recent work (Zhang 2025) gave strong counterexample for the equivalent statement over U 2
f(U) (SU,U) U) and f(U) = IIUT U.U.. 1|
Therefore, even though K-means is widely known to admit nonnegative Murer-Monteiro reformu- lation (Peng & Wei, 2007). there have been only two prior works that actually follow this approach, to the best of our knowledge. Neither of these rigorously guarantee global optimality under Assumption D. The first is the first-order Riemannian method introduced by Carson et al. (2017).
solves the following:
min UEM { (XXT ,UTT 111
(4)
where M' = {U E Raxk : IK, UU UU 1, 1n and U_ max 0} is the (entry- wise) negative part of U, 0 the penalty parameter for U 2 0. Although superficially similar, their approach fundamentally lacks convergence guarantee to a second-orde critical point, due to: (i) their nonsmooth objective; (ii) their use of a smooth penalty, which cannot truly enforce feasibility U 2 (iii) their use of first-order method, which can get trapped saddle point. Moreover, their manifold is geometrically complicated, necessitating an expensive retraction that costs Onn²) time, which prevents their method from scaling to large datasets.
The second work is the nonnegative low-rank (NLR) method Zhuang et al. (024). This is simple projected gradient descent that directly projects U onto the nonnegative spherical constraint and deals the row sum constraint UU 1, = 1, via the augmented Lagrangian method. It a first-order primal-dual method that can only achieve local linear convergence in a neighborhood of its global solution. Like Carson et al. (2017), it is unclear whether there is pathway that this algorithm can LikeoCarson lead to global optimality guarantee, or even to second-order optimality.
Recently, hybrid methods (Wang & Hu 2025 Monteiro et al. 2025 Hou et al.) 2025 Wang et al.) 2023) have been proposed for low-rank SDPs. These approaches handle simple manifold con- straints via Riemannian optimization, while enforcing the remaining constraints through augmented- Lagrangian updates. In contrast, our reformulated K K-means problem requires strict feasibility with respect to both nonnegativity and the simplex-type manifold constraints, since the recovered factor must encode a valid partition, gogmmeted-LLgrangeee or projection based schemes do not preserve this property and would break the structural guarantees which our method relies.
3

---

Published as a conference paper at ICLR 2026
2 BACKGROUND
2.1 SDP RELAXATION OF K-MEANS
Despite the worst-case NP-hardness of the K-means clustering optimization problem 1 common
practical heuristics and relaxed formulations like Lloyd's algorithm (Lloyd] 1982). spectral cluster-
ing (Ng al. 2001 von Luxburg 2007), nonnegative matrix factorization (NMF) (He al. 2011
Kuang et al. 2015; Wang & Zhang, 2012) and SDPs (Peng & Wei 2007 Royer 2017 Fei & Chen
2018 Chen Yang a02 work surprisingly well solving for weal-world data. To explain this
discrepancy between theory and practice, suppose that the data X1. XI...., Xn E R" are generated from standard Gaussian mixture model (GMM)
X; = A4 E, E 1id. O(0,02 for i G
where G denotes the ground truth clusters. Chen & Yang (202b)) proved that the SDP (2) of Peng & Wei (2077) (as - 0+) achieves a sharp phase transition on the separation of centroids for the clustering problem, in any dimension d and sample size n. Let
θº == 40
(
14
Kd n log T3
log n,
(6)
and mmin == minggjies 11 /Axll be the minimum centroid separation. Assume that n/ K is
an integer without loss of generality and consider any a > 0. As soon as the exact recovery becomes possible in the regime min 2 a1)a the SDP approach (2) solves the K K-means problem without clustering error with high probability. For precise statements on the information-thhereee
threshold, please refer to Theorem 3 in Appendix C As immediate consequence of the global optimality guarantee of the K K-means SDP in (2). we deduce that the global solution of the nonconvex low-rank SDP in 3 solves the K-means clustering problem in ① in the exact recovery regime.
Next, from the membership matrix Z, we would like to convert to the cluster label.
Lemma 1. Let Z z E k be the symmetric block-diagonal matrix defined by Zi, 1/|Gal
if E Gk, and Zij = 0 otherwise. Then for any integer T E [K, nl, there is a unique (up to column permutation) U E R such that 7 = UUT, Moreover U can be recovered from any Ú E TRa satisfying Z ÚŪT in n poly(r) time.
For each block-diagonal membership matrix Z, the unique U RT*K in Lemma | iss the associated group assignment matrix, i.e. the k-th column of provides a one-hot encoding of membership in the k-tth cluster.
2.2 CRITICAL POINTS IN CONSTRAINED OPTIMIZATION
The problems considered in this paper are instances of the following
min ((U), M = +U E Raxx AAUUU) + B(U) = =c} UEM
where the linear operators A: R= + Rm and B: Raxr + RM and right-hand hand side c E R together are assumed to satisfy the linear independence constraint qualification (LICQ)
2[ (v)JU B (y) = y=0 0 YU E M.
(8)
In this context, U E Raxr is said to be feasible if it satisfies U E M. The feasible point U is an c-first-order critical point if it satisfies
V 2|A (V)JU +B BBT()
exists
y E RM
s.t.
(9)
and an e-second-order critical point it additionally satisfies
< (V)( A ()).000)) 2 VElLO112 VÙ E T
(10)
over the tangent space of M at the point U, given by Tu M (Ü E Rax A(UÚT + OUT)= +
B(Ū)- = 0}. Under LICQ (8) every local minimum (and hence the global minimum) is guaranteed to be e-second-order critical point (for any 2 0). Unfortunately. there is general-por -urpoose algorithm that is guaranteed to converge to a critical point, due to the need to achieve and maintain feasibility across all iterates.
4

---

Published as a conference paper at ICLR 2026
2.3 SECOND-ORDER RIEMANNIAN OPTIMIZATION
Riemannian algorithms are special algorithms that maintain feasible iterates through a problem- specific retraction operator, and hence able to rigorously guarantee convergence critical points. The basic idea is to improve a feasible iterate U E M by tracing a smooth curve on the feasible set 7: (0,e) - M that begins at 7(0) U and proceeds in direction of descent r(0) = U E To M. In analogy with unconstrained algorithms, a good choice of E Ty M is found through local Taylor expansion
f(((t) f(U) t(grad f(U) + ta (Hess f(U) O(t))
(11)
where grad f and Hess f are respectively the Riemannian gradient and Riemannian Hessian of on the manifold M. Afterwards, we trace the curve y(t) = (Rt(ttt)) using second-order retraction operator Ry: Ty M - M satisfying
d
=Ü, d,
(0) U,
Ry RO() di
(1( dr
T To M,
1=
t=
for all U E M and all Ú E To After choosing step-size 1 so that Unew = 7(1) makes a sufficient improvement over U, repeat the algorithm until it reaches e-second-order critical point satisfying ||grad [(U) S e and Amin (Hess f(U)) 2 -VG, which incidentally corresponds exactly to (9) and [10. Proofs for the following convergence result can be found in Zhang & Zhang (2018); Boumal et al. (2019); AAgrral et al. (202T); we have chosen the simplest but most restrictive settings to ease the exposition.
Theorem 1 (Riemannian cubie-regularized Newton). Suppose that minueM f(U) -o, and that the pullback f f o Ry has Lipschitz continuous Hessian for all U -M. Then, there exists sufficiently large regularizer L such that Uk+1 RUg (Ux) where
ÜR arg min f(U) + (grad f(U). Ū) + (Hess S(U))OT, + - 101 UETU M
converges to an e-second order critical point in O Oce33/2 iterations, independent of dimension.
Each iteration of Riemannian cubic-regularized Newton solves an expensive Newton subproblem. Although it converges in far fewer iterations compared to gradient methods, it is practically competi- tive only when the added cost of solving the Newton subproblem can be offset by the corresponding reduction in iteration count.
3 FORMULATION AND SOLUTION OF K-MEANS AS MANIFOLD OPTIMIZATION
We now explain how we solve (3) using a Riemannian optimization approach. As a first attempt, we can indeed verify that the the constraint set in (3), written
M= {U Raxx, LUTT,, 1,, 1n, tT(UUT)= K}.
(12)
is manifold by checking that (8) holds (cf. Lemma 3 in the appendix). In fact, directly applying Riemannian optimization techniques results K- K-means algorithm very similar to the one proposed
in Carson et al.] (2017). The immediate and critical difficulty with this approach is the lack of an efficient retraction operator, which must be called at every iteration to keep iterates feasible U M. For example, Carson et al. 2017)) used complicated exponential retraction that costs O(n²) time, hence bottlenecking the entire algorithm and preventing it from scaling to large n.
Instead, our first contribution in this paper is to reformulate 3 by establishing submersion from the product manifold M x Orth(r) to M, where
V {v Rax((-)) 1, 0 tr(VV = -K-1] Orth(r) {Q xx QQ = 1r}.
In words, V is projected hypersphere and Orth(r) is the of x orthonormal matrices.
5

---

Published as a conference paper at ICLR 2026
Let (pVV Q) VQ. where V == [in V]with i, == (1/vn)ln Then M ((M). Moreover, the Jacobian Dp: T M + T M is surjective for all (V, Q) E M, i.e., is submersion.
Theorem 2.
Having established the submersion property of ip it standard result that every c-second order point of M is also an ce-second order point on M for some constant rescaling factor c; see e.g. Example 3. 14 and the surrounding text in Levin et al. (025). Therefore, to solve (3), equivalently solve
min (C, (C, a Σ log (o)((P)(V (VIQ)EM
(13)
where C XX is the (negative) data Gram matrix, and Piy the (1.j)-th element of the operator in Theorem 2 A basic but critical benefit of the reformulation 13 manifold M admits simple second-order retraction via its Euclidean projection (Boumal 2023
is that the product
Sec. 5.12)
Riv.o) Q) (Vojyv(++ + +V) Projonn()(Q++0))
where
V - n n V. IV - -
Projy(V) VK
and Projonth(r) (Q) = (QQ
It is easy to check that the retraction above costs just O(nr + r²) time to evaluate. In Section E.2 we give explicit expressions for the Riemannian gradient and Hessian and explain how they can be computed in O(rr + r²) time.
The appearance of the logarithmic penalty in (13) presents two difficulties. First, as a practical concern, any algorithm for 13) must begin at a strictly feasible point (Vo. Qn) E . that additionally satisfies yo (V,, Qo) > O. In Section E.4 we provide a good strictly feasible initial point, and prove that points exist only if the search rank is over-parameterized as > K. Second, some special care needed to rigorously apply the guarantees from Section 2.3 given that the penalty <VVV Q) is Lipschitz only when restricted to a closed and strictly feasible subset; see Section E.5 for details.
Together, these ingredients allow to apply Riemannian gradient descent (Boumal al. 2019) to (13 to compute an e- e-dirst-order critical point in (n/c) poly(r, d, K) time. In practice, the algorithm often converges an e-second-order critical point, though this not rigorously guaranteed without carefully-tuned noise perturbation. Alternatively, we can apply the seajjugttegrrdieeni (CG) variant of the Riemannian trust-region algorithm (RTR), a general-purposs solver available in packages like MANOPT (Boumal et al.) 2014) or PYMANOPT (Townsend et al., 2016). to guarantee convergence to an e-second-order critical point. Unfortunately, in our experiments, we observed that all of these algorithms experience unsatisfactorily slow convergence, due the severe ill-conditioning introduced by the logarithmic penalty.
Instead, our best numerical results were obtained by the Riemannian cubic-regularized Newton (Theorem l.. Our key insight is that the algorithm can be implemented with just O(er²) time per-iteration, by exploiting the underlying Mlec--diggonll-luusloo----- Hessian. To explain, our core difficulty is to efficiently solve the Newton subproblem
structure of the Riemannian
min g ppp 2 Hp |Pll
(14)
Ap=0
where g and H denote the vectorized Riemannian gradient and Hessian respectively, and A imple- ments the tangent space constraint (V, Q) E T(v:) M. We can verify that the subproblem contains m(r 1) + 7² ((r) variables and is subject m + r(r 1)/2 O(r²)) constraints. Given that the subproblem has only linear constraints, its local minima must always satisfy the first- and second-optimality conditions ( and (T).. which read
H A1 A
A
ppl
- A= Ipl.. (T AI) 2 0 for all & satisfying Ag = 0. 1Pl,,
The following result be viewed Riemannian extension of known results on the approximate minimization of regulaaegdl subproblems. It shows that, with sufficient regularization L, the global minimum corresponds to the unique second-order critical point.

---

Published as a conference paper at ICLR 2026
Lemma 2. Let A have full row-rank (i.e. AAT -Amiu the parameterized solution
0) and let Xmin minjel 11, HÇ. For
D 0
A A
p(A) =
A
0
well-defined and [|x( is monotonously decreasing with respect to X.
The same lemma also suggests solving the Newton subproblem by simple bisection search, ef. (Cartis) et al. 2011 Sec. 6.1). Indeed, the solution is just p(Aop)). where Acpt is the solution to the monotone
equation 2A- [|)] (via Lemma 2). Thus, we pick very small Ab R -Aain such that ||p( (n(Ab)] > 2/mh I., very large Auh such that 2Auh/T. >> [|p[ II. and then perform bisection until 27ugI T. (1Ipp()op))] is approximately found. For each X, if 2 T.]p(A)]],, then we increase A;
otherwise, we decrease A.
The main cost of the bisection search is the computation of p(A), which naively costs O(n"r")) time. For our specific problem, we explain in Appendix [ how a block-diagonal blck-ddaggnal-plus-oowwzza the Hessian H reduces the computation cost just n poly(r, d) time. Applying Theorem I shows that the overall method computes an e-second-order critical point in (n/e!.)) pplly(r, d, K) time.
structure in
4 NUMERICAL RESULTS
In this section, we showease the superior performance of our proposed Riemannian second-order method for clustering on both synthetic Gaussian mixture models (GMM) and real-world mass cytometry (CyTOF) datasets. Compared to existing oftte-of-thhaatt methods, such as the nonnegative low-rank (NLR) factorization (Zhuang al., 2024) and prior Riemannian K-means algorithms (Carson et al. 2017), our approach achieves faster convergence, higher clustering accuracy, and more reliable recovery of ground-truth cluster memberships. These results highlight the convergence and accuracy advantages of second-order methods when they can be implemented with per-iteration costs of just O(n) time. The implementation details are deferred Appendix H
Datasets. We conducted experiments on both synthetic and real datasets. The synthetic data was generated from standard K componnen, d-dimensional Gaussian mixture model (GMM), with centroids placed at simplex vertices such that their separation equals - where θ is the information- theoretic threshold for exact recovery in (6), and y controls separation. The real dataset came from mass cytometry (CyTOF) (I (Levine et al. 20155 Weber 2015). It consists of 265,627 cell protein expression profiles across 32 markers, labeled into 14 gated cell populations. Following /Zhuang et al.. (2024), we uniformly sample 1,800 cells from K 4 unbalanced clusters (labels 2, 7, 8, 9) from individual 1 for our experiment.
Global optimality at second-order critical points (validation of Assumption D. Figure 1 shows the convergence behaviors of loss function (13) for GMMs (n = 500, 7 = 1.2, u 0.0T) with 50 randomized initializations. We consistently observe that: (i) the loss value steadily decreases over iterations and converges rapidly near the globally optimal point; (ii) the Riemannian gradient norm dynamics suggest that algorithm initially attempts escape saddle points (with increased gradient norm) and eventually converges to second-order local optimality, where zero-loss is achieved, indicating global optimality. To verify second-order local optimality, we also plot the minimum eigenvalue of the Riemannian Hessian. This provides strong numerical evidence that near-second- order critical points are rear-globally optimal, as posited by Assumption
Benchmark on real world data. Prior studies on mass cytometry (CyTOF) and computer vision (CIFAR-10) datasets identified the nonnegative low-rank (NLR) factorization (Zhuang et al.) 2024) as the most reliable clustering solver, attaining the lowest average mis- mis-clustering error and the tightest variance compared to classical baselines such spectral clustering (SC), nonnegative matrix factorization (NMF), and K-means++ (Arthur & Vassilvitskii 2007) (KM++). Our algorithm optimizes the same nonnegative low-rank model, so it inherits this reliability. Because it applies second-order Hessian updates rather than first-order gradients, it refines each iterate more thoroughly and therefore recovers the ground-truit membership matrix more accurately. Figure 3 illluttrates this on CyTFF: both methods keep mis- mis-clustering near zero, yet our solver achieves a smaller Frobenius gap to the oracle solution. The experiment was repeated 50 times on random subsamples of size 1800. An additional experiment on the CIFAR10 dataset can be found in Appendix G
7

---

Published as a conference paper at ICLR 2026
Comparison with NLR. Next, we compare method directly to the nonnegative low-rank (NLR) factorization. Figure 2ishows experimental results for GMM with n = 100, ? = 0.8, μr = 0.1, and with varying Main observations: (i) Each Newton step is solved in (1)) time, matching the theory in Section 3 (ii) A Newton step is about 30-100 times costlier than a single NLR update. This is to be expected because the Newton step solves several linear systems, while NLR performs only single matrix-vector product. (iii) Our solver reaches the optimum in hundreds of iterations, whereas NLR needs tens of thousands. The ordens-of-mannitude reduction in iterations more than offsets the costlier step, so wall-clock time drops by factor of two to four.
Comparison with prior Riemannian K -means methods. We evaluate the clustering performance of the algorithm proposed by Carson et al. (2017) to solve the penalized formulation (4) on the CyTOF data. Figure 4 presents the performance of this first-order manifold method. Unfortunately, we were unable to identify sequence of À that would produce acceptable clustering results. While increasing the penalty parameter improves feasibility, also degrades clustering performance. This highlights the difficulty in solving (4) using the method of Carson et al. (2017), as it struggles to balance strict constraint satisfaction with objective optimality in the -meass problem.
Comparison with classical Riemannian algorithms. As discussed in Section Problem 13 can be solved with CG or RTR as the gradient, Hessian, and retraction are all available. Nevertheless, both CG and RTR perform poorly because the log- oog-barieer induces extremely ill ill-conditooned landscape. To illustrate this, we solve (33) on GMM data using PYMANOPT's implementation of CG and RTR. For unable to tune the method to produce meaningful solution, its updates frequently lead to infeasible points. While RTR can converge to a solution comparable to ours, it requires significantly more iterations and time, shown in Figure
CG,
Effect of hyperparameters. Our method displays a sharp phase transition with respect to the regularization parameter p. Below critical threshold, the algorithm consistently converges to
Global Suboptimality
Gradient Norm
Minimum Hessian Eigenvalue
10a
10
10a
10
10-6
10-1
10
250
Iterations
Iterations
Iterations
Figure 1: Local convergence to second-order critical points yields global optimality. In the GMM setting, where ground-truth partitions can be planted, we consistently observe local convergence to global optimum, yielding zero clustering error. This provides strong numerical evidence that rear-second-order critical points are near-globally optimal, as hypothesized in Assumption |
Mis-clustering Error
Mis-clustering Error
Per-iteration Time Cost (sec)
Ours
10P
Ours
0.6
0.6
NLR
NLR
0.4
0.4
10-2
10-3
0.2
Ours
0.2
10l
10a
10a
10
20
10a
10
Iterations
Time Elapsed (soc)
Sample Size
NLR on GMM. Our second-order method reaches optimality in 152 iterations, while NLR needs 80k. Even though each second-order iteration costs x5 25-100 NLR steps, the total runtime is still two to four times shorter. (Left and middle) clustering accuracy log iterations and linear time. (Right) per-iteration time sample size
Figure 2: Comparison with previous
8

---

Published as a conference paper at ICLR 2026
Mis-clustering Error
Frobenius Distance to Z*
8
0.4
1.2
0.3
1.0
0.8
0.2
g
0.6
0.1
0.4
0.0
Ours
NLR
SC
NMF
KM++
Ours
NLR
Figure 3: Real-world benchmark on CyTOF data. We compared our method to NLR, the previous tfae-offthhe-ar, as well as classical benchmarks SC, NMF, and K M+++ Our method and NLR achieve the most consistently accurate clustering, with the smallest variance and the fewest outliers (left), but we outperform NLR in ground truth recovery (right).
Mis-clustering Error
Frobenius Distance 2"
Norm of the Negative Part
2.0
1.5
0.4
1.5
0.2
0.5
0.0
Ours
Ag A4 Ours
7
Ag
Ag
24
A Ag A4
Figure Comparison with prior Riemannian K-means method Carson et al. (2017) on real-world data. Each run is warm-started from the previous and the penalty is stepped through A = 0, _0a 105, 107. However: (Left) average mis-clustering exceeds 30%; (Middle) the recovery error 2 - 1= remains large; (Right) the infeasibility ||U_|| never vanishes. Our Riemannian method, shown for reference, enforces U__ = by design and achieves near-zero error in both metrics.
Global Suboptimality
Gradient Norm
106
10'
10
10
109
10
10-3
10-8
Ours
10-11
10
RTR
1000
1500
2000
5000
10000
15000
20000
Figure 5: Comparison with classical Riemannian Trust Region (RTR) on GMM. Our method drives both loss and gradient norm machine precision around 360 iterations. In contrast, RTR stagnates for over 2 1k iterations due to the extreme ill-conditioning induced by the log penalty.

---

Published as a conference paper at ICLR 2026
optimal solutions and remains robust variations in μ.. However, once ja exceeds this threshold, the method fails to yield meaningful results. This behavior is illustrated in Figure 6 where we used GMM synthetic data with four clusters and separation y 0.8. Figure also shows the clustering errors evaluated across different search ranks r. Clustering performance is largely insensitive to when js is small, the smallest feasible K 1 is recommended in practice. However, when is large enough to trigger the the phase transition, increasing r can help the algorithm escape spurious regions. A brief discussion the choice of hyperparameters provided in Appendix
Mis-clustering
Mis-clusten Mis-clustering Error
Mis clustersgg Error
0.05
0.6
0.04
0.03
0.02
0.5
7
5
Figure 6: Dependence on hyperparameters. (Left) Error exhibits phase transition as penalty parameter gets too large. (Middle): Errors are insensitive to the search rank when is chosen appropriately. (Right): When too large and the algorithm becomes trapped in local minima, increasing can leads to significantly better results.
Robustness to mis-specified cluster number. Figure 7 shows the performance of our method under mis-specified cluster numbers in the GMM setting with K 4, n 1000, and 7 0.8; The mis-specified cluster number is set to be Kmis = 3 and Kmis 5. In both setups, we observed that our method exhibits strong robustness in statistical accuracy for label recovery. We used normalized mutual information (NMI) to quantify the agreement between true and recovered labels, as the two sets labels have different categories. Moreover, we observed that our method locally converges to the global optima in similar manner as the K = 4 case (cf. Figure
Normalized Mutual Information
Global Suboptimulity
Gradient Norm
0.85
10a
10a
10-1
10a
0.80
10
10
10-7
0.75
10
10
Ours
NLR
SC
XM++
100
200
100
200
1.00
0.95
10a
10a
10-1
0.90
10a
10
0.85
10-1
10
0.80
10-10
10
Ours
SC
XM-++
500
7500
Iterations
Iterations
Figure 7: Robustness to missspeciiication. Clustering performance and convergence behavior when the number of clusters is under-estimated (top) and over-estimated (bottom)
10

---

Published as a conference paper at ICLR 2026
ACKNOWLEDGMENTS
X. Chen was partially supported by NSF DMS-2413400 2413404 and a gift from the Simons Foundation. R. Zhang was partially supported by NSF CAREER Award ECCS-2047462 and ONR Award N00014- 24-1-2671.
REPRODUCIBILITY STATEMENT
The experimental setups are described in detail in Section and Appendix G does not encompass every experiment reported in the paper, all of them can be readily reproduced based on the descriptions.
While the provided code
LLM USAGE DISCLOSURE
A large language model (LLL)) was used solely to polish the writing and improve clarity of expression. No part of the research ideation, discovery, or substantive content was generated by the LLM.
REFERENCES
P-A Absil, Jochen Trumpl, Robert Mahony, and Ben Andrews. All roads lead to newton: Feasible second-order methods for constriiee optimization, 2009. Technical Report.
Naman Agarwal, Nicolas Boumal, Brian Bullins, and Coralia Cartis. Adaptive regularization with cubics on manifolds. Mathematical Programming. 188:85-134, 2021.
Daniel Aloise, Amit Deshpande, Pierre Hansen, and Preyas Popat. NP-hardness of Euclidean sum-of-squares clustering. Machine learning, 75(2))245-248, 2009.
David Arthur and Sergei Vassilvitskii. k-means++; the advantages of careful seeding. In Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA '07, pp. 1027- 1035, USA, 2007. Society for Industrial and Applied Mathematics. ISBN 9780898716245.
Afonso S Bandeira, Nicolas Boumal, and Vladislav Voroninski. On the low-rank approach for semidefinite programs arising in synchronization and community detection. In Conference on learning theory, pp. 361-382. PMLR, 2016.
Srinadh Bhojanapalli, Behnam Neyshabur, and Nati Srebro. Global optimality of local search for low rank matrix recovery. Advances in Neural Information Processing Systems, 29, 2016.
N. Boumal, B. Mishra, P.-. Absil, and R. Sepulchre. Manopt, Matlab toolbox for optimization on manifolds. Journal of Machine Learning Research, 15545):445511459, 2014. URL https wwww.manopt.oo org
Nicolas Boumal. An Introduction Optimization Smooth Manifolds. Cambridge University Press, 2023. doi: 10.1017/9788 1009166616 166164 URL /ttttps://ww.. URL/https:////
nettbbook
Nicolas Boumal, Vlad Voroninski, and Afonso Bandeira. The non-convex burer-monteiro approach works on smooth semidefinite programs. Advances Neural Information Processing Systems, 29, 2016.
Nicolas Boumal, Fierre-Antoine Absil, and Coralia Cartis. Global rates of convergence for nonconvex optimization on manifolds. IMA Journal of Numerical Analysis, 39(1):1-33, 2019.
Nicolas Boumal, Vladislav Voroninski, and Afonso S Bandeira. Deterministic Guarantees for Burer- Monteiro Factorizations Smooth Semidefinite Programs. Communications on Pure and Applied Mathematics, 73(3)5581-608 2020.
S.P. Boyd and L. Vandenberghe. Convex Optinization. Number pt. 1 in Berichte über verteilte messysteme. Cambridge University Press, 2004, ISBN 9780521833783, URL lttps btttpp::/bbolke bookss google cow./bookssiiddmmnob.b.....

---

Published as a conference paper at ICLR 2026
Samuel Burer and Renato D. C. Monteiro. A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization. Mathematical Programming, 95(2):329 -577 2003. ISSN 1436-4646. doi: 10.10077/110107002--55333 URL https://doi. org/ /t.g//0..00 1007 610107-002-0352-3
Richard H Byrd, Jean Charles Gilbert, and Jorge Nocedal. A trust region method based on interior point techniques for nonlinear programming. Mathematical programming. 89:149-185, 2000.
Richard H Byrd, Jorge Nocedal, and Richard A Waltz. Knitro: An integrated package for nonlinear optimization. Large-scale nonlinear optimization, pp. 35-59, 2006.
Timothy Carson, Dustin G. Mixon, Soledad Villar, and Rachel Ward. Manifold optimization for k-means clustering. In 2017 International Conference on Sampling Theory and Applications 200.://AAMMA.21...822222
(SampTA), pp. 73 73-77, 2017. doi:
Coralia Cartis, Nicholas I. M. Gould, and Philippe L. Toint. Adaptive cubic regularisation methods for unconstrained optimization. Part I: motivation, convergence and numerical results. Mathematical Programming, 127(2):245-22) 295, 4 2011. ISSN 1436-4646. doi: 00.1007/110107-000-222222
Xiaohui Chen and Yun Yang. Hanson- Wright inequality in Hilbert spaces with application to K-means clustering for -on-Euclidean data. Bernoulli, 27(1))586 614, 2021a. doi: 10.3150/ 20- BEJ1251.
Xiaohui Chen and Yun Yang. Cutoff for exact recovery of Gaussian mixture models. IEEE Transac- fions on Information Theory, 67(6))4223 4238, 2021b. doi: 10.1 00.10099TT2222233333555
Sanjoy Dasgupta. The hardness of k-means clustering. Technical Report CS2007-0890, University of California, San Diego, 2007.
Peter JC Dickinson and Luuk Gijben. On the computational complexity of membership problems for the completely positive cone and its dual. Computational optimization and applications, 57: 403-415, 2014.
Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, Jakob Uszkoreit, and Neil Houlsby. An image is worth 16x 16 words: Transformers for image recognition at scale. In International Conference on Learning Representations, 2021. URL, http:ttttss///erreeet nett/Eorum?id--icbbatt
Yingjie Fei and Yudong Chen. Hidden integrality of sdp relaxations for sub-gaussian mixture models. In Sébastien Bubeck, Vianney Perchet, and Philippe Rigollet (eds.), Proceedings of the 31st Conference On Learning Theory, volume 75 of Proceedings of Machine Learning Research, pp. 1931-1966 PMLR, 06-09 Jul 2018. URL httpttts:// html
mlr ppressv75/ffilll
Rong Ge, Jason D Lec, and Tengyu Ma. Matrix completion has no spurious local minimum. Advances in neural information processing systems, 29, 2016.
Rong Ge, Chi Jin, and Yi Zheng. No spurious local minima in nonconvex low rank problems: A unified geometric analysis. In International Conference on Machine Learning, pp. 1233-1242. PMLR, 2017.
Zhaoshui He, Shengli Xie, Rafal Zdunek, Guoxu Zhou, and Andrzej Cichocki. Symmetric non- negative matrix factorization: Algorithms and applications to probabilistic clustering. IEEE Transactions on Neural Networks, 22(12):117-2133, 2011.
doi: 10. 20.11999NNN2011 20122.212..2..
Di Hou, Tianyun Tang, and Kim-Chuar Toh. A low-rank augmented lagrangian method for doubly nonnegative relaxations of mixed- binary quadratic programs. Operations Research, 10 2025. ISSN 0030-364X. doi: 10.1 1287/opre.2024.1 137.
Wen Huang, Meng Wei, Kyle A. Gallivan, and Paul Van Dooren. A riemannian optimization approach to clustering problems. Journal of Scientific Computing, 103(1):8, Feb 2025. ISSN 1573-7691. doi: 10.10077/10915-005--0000
12

---

Published as a conference paper at ICLR 2026
V. Kalofolias and E. Gallopoulos. Computing symmetric nonnegative rank factorizations. Lin- ear Algebra and its Applications. 436(2):421-435 2012. ISSN 0024-3795. doi: 10.
(016/).la.2011103.011 URL https://www sciencedirect. httm///cience//ticcleloo/oo
$0024379511002199 Special Issue devoted to the Applied Linear Algebra Conference (Novi Sad 2010).
Etienne Klerk. Aspects of Semidefinite Programming. Springer New York, NY, 2006. ISBN 978-0-306-44819-2. doi: 10.1007/b105286.
Da Kuang, Sangwoon Yun, and Haesun Park. SymNMFF: nonnegative low-rank approximation of similarity matrix for graph clustering. Journal of Global Optimization, 62:545-574 2015.
Eitan Levin, Joe Kileel, and Nicolas Boumal. The effect of smooth parametrizztions on nonconvex
optimization landscapes. Mathematical Programming, 209(1))63-111, 01 2025. ISSN 1436-4646. doi: 10000///000700024000555555
Jacob Levine, Erin F. Simonds, Sean Bendall, Kara L. Davis,
D. Amir, Michelle D. Tadmor, Oren Litvin, Harris G. Fienberg, Astraea Jager, Eli R. Zunder, Rachel Finck, Amanda L. Gedman, Ina Radtke, James R. Downing, Dana Pe'er, and Garry P. Nolan. Data-driven phenotypic dissection of aml reveals progenitor-like cells that correlate with prognosis. Cell, 162(1):184-11)))) 2015. ISSN 0092-8674. doi: 10. 10..00016f..el2221555....4 URLhttpps h/tp://wwwwwccccc www...ccieceeirectt
com com/ httm/scienceeeee http:///eeaaccacaaeea//////2//2/
Stuart Lloyd. Least squares quantization in PCM. IEEE Transactions on Information Theory, 28:
129-137, 1982.
Jan R. Magnus and Heinz Neudecker. Matrix Differential Calculus with Applications in Statistics and Econometrics. Wiley Series in Probability and Statistics. John Wiley, third edition, 2019. ISBN 9781119544 202. doi: 10.1002/9781 19541219.
Bamdev Mishra and Rodolphe Sepulchre. Riemannian preconditioning. SIAM Journal on Optimiza- lion, 26())66356600, 2016.
Renato D.C. Monteiro, Arnesh Sujanani, and Diego Cifuentes. A low-rank augmented lagrangian method for large-scale semidefinite programming based on a hybrid convex-nonconvex approach, 2025. arXiv:2401.12490.
Katta G Murty and Santosh Kabadi. Some NP-complete problems in quadratic and nonlinear programming. Mathematical Programming, 39:117-129, 1987.
Andrew Ng, Michael I. Jordan, and Yair Weiss. On Spectral Clustering: Analysis and an algorithm.
In Advances in Neural Information Processing Systems, pp. 849-856. MIT Press, 2001.
Jorge Nocedal and Stephen J. Wright. Numerical Optimization. Springer Series in Operations Research and Financial Engineering. Springer New York, NY, 2nd edition, 2006. ISBN 978-0-387- 40065-5. doi: (0007/978-0-337-4006666
Jiming Peng and Yu Wei. Approximating K K-means-type clustering via semidefinite programming. SIAM J. OPTIM, 18(1):186-205, 205, 2007.
Wei Qian, Yuqian Zhang, and Yudong Chen. Structures of spurious local minima in k-means. IEEE Transactions on Information Theory, 68(1):395-422, 2022. doi: 10.100991T.2221131224444
Martin Royer. Adaptive clustering through semidefinite programming. In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (cds.), Advances in Neural Information Processing Systems 30. pp. 1795-1803. Curran Associates, Inc., 2017.
James Townsend, Niklas Koep, and Sebastian Weichwald. Pymanopt: A Python Toolbox for Opti- mization on Manifolds using Automatic Differentiation. Journal of Machine Learning Research, 17(137):1-5, 2016. URL htt.httt:://// org/ wwtpg/ppper///////---- 177 html
Ulikee von Luxburg. A tutorial on spectral clustering. Statistics and Computing, 17(4):395 416. 2007.
13

---

Published as a conference paper at ICLR 2026
Jie Wang and Liangbing Hu. Solving low-rank semidefinite programs via manifold optimiza- tion. Journal of Scientific Computing, 104(1):33, 2025. ISSN 1573-7691. doi: 10.1007/
s10915-025-02952-8
Yifei Wang, Kangkang Deng, Haoyang Liu, and Zaiwen Wen. decomposition augmented lagrangian method for low-rank semidefinite programming. SIAM Journal on Oplimization, 33(3):1361-1390 2023. doi: 10.1137/22M14745599
Yu-Xiong Wang and Yu-Jin Zhang. Nonnegative matrix factorization: comprehensive review.
IEEE Transactions knowledge and data engineering, 25(6):133--1355 2012.
Lukas Weber. Clustering benchmark data: 32-dimensional data set from Levine al., 2015. URL GitHub Repository.
Junyu Zhang and Shuzhong Zhang. A cubic regularized newton's method over riemannian manifolds,
2018. AXXi:1800.005560
Richard Y Zhang. Nonnegative low-rank matrix recovery can have spurious local minima, 2025. arXiv:2505.03717.
Yubo Zhuang, Xiaohui Chen, and Yun Yang. Wasserstein K-means for clustering probability distributions. In Advances Neural Information Processing Systems, 2022.
Yubo Zhuang, Xiaohui Chen, Yun Yang, and Richard Y. Zhang. Statistically Optimal K K-means Clustering via Nonnegative Low-rank Semidefinite Programming. In The Twelfth International Conference on Learning Representations, 2024.
14

---

Published as a conference paper at ICLR 2026
A RELATIONSHIP BETWEEN SDP FORMULATIONS OF K-MEANS CLUSTERING: STANDARD AND INTERIOR POINT VERSIONS
The standard K-means clustering for data in Rd has two formulations in the literature: (1) centroid- based optimization
min P1 ERI
min kE(K]
|X - 8 12:
and (ii) partition-based optimization
ΣΣ |X Xx112.
min
U GL=(n] kE1
where Xx = 1Ga= G GG X, is the empirical centroid of cluster Gs. Formulations (i) and (ii) are known to be equivalent, cf. Zhuang et al. (2022) Eqn. (1)) or Qian et al. 2022 Appx. A). Using the
parallelogram law in Zhuang al. (2022 Eqn. (5))
Σ
| x,l12 2|
Σ |X Xx
KJEC
LEGx
may write the partition -based objective function
min
Σ LJEG2
-
U Gy [+]
2
Next, expanding the pairwise squared Euclidean distance and dropping Σi= depending on any partition GI...., GR), we arrive at (T). in the form of maximizing the total intra-cluster similarity in terms of the Gram matrix {(x,, X
(no longer |X,13
For general data without a likelihood derivation as in Chen Yang 2022b). we can replace the Re-inner product with any (positive semidefinite) kernel k: X' X + R and consider the ker- nelized version of K K-means clustering that involves data XI...., Xn only via their Gram matrix {k(X,, j= Thus, our manifold formulation of this paper carries over to the general kernel method setting with possibly nonlinear boundary structure.
Next, convexify the K K-means problem (I) into SDP. Note that each partition GI...., GG via one-hot encoding is equivalent to an assignment matrix Ha*K (up to cluster relabel) where the latter is binary matrix with exactly one non-zero entry in each row, i.e. H₁ 1 if E Gk. With this reparameterization, one can write [T)
as a mixed zero-one integer program:
(XXX HBH H1 = 1,)
max HE(0,I)axK
Now, applying the change of variables Znxn HBH and noting that the membership matrix Z and assignment matrix H are not one-to-one, we relax the K-means problem by preserving the key properties of Z as the following constraints:
z O,
tr(Z)= K.
Z1, In-
Z 0.
which no longer depend on the assignment matrix Then, we arrive at the standard SDP relaxation for K- K-means clustering, cf. Peng & Wei [2077 Equ. (13)) or Chen & Yang 2021b Eqn. (11)):
ZER max { (XXX] Z) Z 0, tr(Z) K, Z1, 1,, 2 2220 0}.
(15)
In practice, the elementwise nonnegativity constraint Z 2 is almost always enforced by logarith- mic barrier. This means that we can make (15) more explicitly in the form
max Z
(XXT 2) +
log (og(Z,,)+ 2 0, tr(2) 0, K, 21, = 1,
(16)
+11=1
which is precisely how any practical interior-point solver would solve the original SDP in 15 The barrier cost is the actual objective used internally, with y set to reflect the solver's target accuracy, cf.
Boyd & Vandenberghe (2004) Chapter 11) or Nocedal & Wright (2006) Chapters 14 & 19). In other words, we take the standard SDP in (15). and make explicit the logarithmic penalty that is already implicit in how such an SDP is actually solved.
15

---

Published as a conference paper at ICLR 2026
B
TIGHTNESS OF PROPOSED FORMULATION
As shown in Carson et al. 2017). the combinatorial K-means problem ① is exactly equivalent to {XXX,, , U 0,00 1, 1,, = -
max UER"XR KR
Formulation replaces the hard constraint U > 0 with log barrier, and then relax U Ix into Lr(UUT) K. Compare to 2.. the containment (UU U 0) c
0, 2 0) is
strict in general, replacing UC by Z yields relaxation. Chen & Yang (002b)) showed that
the resulting relaxation from to (2) is exactly tight above the statistical separation threshold: its
optimizer factors as Z UU with U 2 0, and the sparsity pattern of U recovers the true cluster labels. Since
is tighter than it must also be tight in this regime.
C INTORMATION-THEORTE THRESHOLD FOR EXACT RECOVERY
The following theorem a precise statement of the iooommation-thorettee threshold of Chen & Yang (2021b).
Theorem 3 (Average-case phase transition for exact recovery). Let a > 0 and
mmin
min
||
p
ISIRRKK
be the minimum centroid separation. Suppose that data XI..... ... Xn, are generated from the Gaussian mixture model (5) with equal cluster size |Gil -..-= |G dichotomy.
Then we have the following
1. K 1/ log log(n) and O O (1 G then there exist constants C1, C2
depending only on such that, with probability at least 1 - C1 Clooggn) C2 ,, the SDP 2
μ 0+ (cf. Appendix for the equivalence of two SDP formulations) has unique solution that exactly recovers the true partition G ...
2. If K S log n and Omin I (1 - a)O, then there exists a constant C3 depending only such that
C3K
inf
sup
P(ak GR + G 2
G.....GK ((n,,
where the infimum taken over all possible estimators (G1,... GK) for (Gi Gk)) and the parameter space is defined as
||; - pul θ, Vj, k E [K] + k)
E(n, K, = (GL....., GK. pk)
(1
| (1 5m)m
with 8, = C. MVKK log(n)/n /nn///fr/ffff some large enough constant C > 0.
16

---

Published as a conference paper at ICLR 2026
D PROOFS
Proof of Lemmal 7
Note that the membership matrix Z associated to partition Gi,... Gx contains diagonal principal submatrix of rank K. The lemma follows from Theorem 4 in Kalofolias & Gallopoulos (2012).
Proof of Lemma Let N denote the null space basis of A, such that AN and N Then, we have Amin = Amin (N HN) and p(A) = (A( where (N HN + XD) = Then, À >> -Amin implies that N HN + Al > so p()) (N HN + AT) g always well-defined. Moreover, [|p())]| = [[(()] is monotonously decreasing because all the eigenvalues of N HN AI are strictly positive and increasing with X.
Proof of Lemma 5
To prove the first statement, note that the implication
Rack A AUUT U R 20
= -
= U T Raxk A 1, A |1 = 20 20
straightforward To see the converse, note that U a (doubly) stochastic stochaasti matrix, tr(UUT) K, and ruUk(UUT)) K, thus all the K eigenvalues of UUT are i. U E St(n, K).
For the second statement, it is trivial that
1
K G2 = [n]
G
1GK
A
G1
VIGK
k=1 U RARK A AUU = ,,
The converse is also true. Let us denote Ga == supp(us). Observe that G.n & for all i + j since 2 0 and u. u; 0 for all + j. Then UU 1, = I, implies u.u 1c. la, and L G+ = [n]. Since || = 1, we have that us = lo./V/Gul-
Finally, we will prove the third statement. Let U be group assignment matrix as defined by 35 Note that for all V E To MM n Cv. V must satisfy (UV + VUT),, = O., (U, V) 0. and
Wijj 2 0, Vuu.
Define [n] >> [K] to be the group assigning function, i.e. A(i)
Ek= KIf + 0}. then
1
UV
((,(())
1
(,A(()
22,((1)
A())
ViGAm
v GAmi
v G
1
11,((2
ViGAil
21,A(n)
V/Gamo) Un,)(n)
//Ga(m)
Observe that
(i) A(I) = A()) + i GAG) + GEGA
(ii i $ GAO) + (u,,AG) 0 (),A(() = 0) = (Vi, 2 0A 0)A(0) 20)
(iii) (U.V) =0 - Vj VjEln) [n), Eiegro 2g,(() 0.
17

---

Published as a conference paper at ICLR 2026
Denote (UVT + VUT)LL ), then
wy Σ
Bs.a(1)
GiGA
vE(ni
Σ
1
1
(LA) +
(g,A(()
VIGAGA
VlGA(o)
VEGA()
Vi,A +
(i)(())
VIGAG)I
V(a(())
By (i) and (iii), we know that
Σ
((,00) +
VGa ()A(() = VIGAssllsSGG)
tEGAG)
By (ii), we know that
Σ GGG
(LA()) +
VGG (IA(() 2 0.
VGGA
Thus we can write
Wj V GAGi USAG)) R2 for some Rý 0.
Next, use proof by contradiction. Suppose there exists V E Ty MCC such that V * 0, then there must be both positive and negative entries in V to satisfy V) = O. This implies that there the only entries that can take negative value, To satisfy (U, V) O, there must exist j2 GA(i) such that A22.((5)) > 0. Then for such j2.
exists E [m] such that ,(I)) < 0 since VJ,A())
Wy VIGAGs)lu (/. R2 V/GAGG) A())) + R2, 0.
This contradicts (UVT VUT)I = O. Therefore, Tr M Cu (0).
Lemma 3. If K 2 2, then the set M defined in (22 satisfies linear independence constraint
qualification (LICQ) for all U E Ruxr and is therefore a smooth submanifold of Raxr.
2020 the set M = Raxr, (A, Proof of Lemma 3 Following Boumal et al.
by for all i} is smooth submanifold of Raxo
LICC holds for all U E M,,.e., that En = 0 verify that
if and only y 0. For the definition in (22
Σ yi A.
111y T + y1 T)
I1(1y + y1 |u1 2(yy + y1 UU, v
[ny] UU+ Jy 2( 2 4yo ") + K
2 2[(1 (((1]T))+22oll 7) + vol n||vl + (K 2)
2(1 y 40) + (K- - 2)
2 2(1 y +
So YAA = 0, then from the last line we conclude that 1 1 = and y = 0, so (v. 30) = 0
as claimed. The third line above because ((1y" + y1 ),, = 2 UU 1) 2(y 1) and
[1yy + y1 yyiiaa {
tr((1y + + (1y' +y1 )
=rr[y T "U y1 + y1 T T1"T"" +2- 1y UU T 1y
=u(1y UUT y1 +1 11y +2- 1y 1yy 1y
+ ny y (1
18

---

Published as a conference paper at ICLR 2026
Lemma 4 (Interior point construction for M,). Given K N. for any such that K, for large enough n, we have the following two cases:
Case 1: n 0 (mod r)
Denote 9 n/r, let Up (x --y) + y11 where
1
1 (r
y -vr- (Vr-I-VK-) VK -1
Then = -1/v & Ug is an interior point of M,
Case 2: n # 0 (mod r) Let us denote q = [n/r] and p = n mod r. Construct the block matrix B E Rrxr;
B B1,2
B =
B2,1
B2
where
B1, (2 y) y1p17
B1, 11,p11pp
B2 (w - 2) 21, 21r-p1,
B2.1 -yrr-p1p
The coefficients z, y, and w depends K and r. They will be specified in the proof. Then
U == Ln 01, B)
an interior point of M,.
Proof of Lemma For general large enough n, is either divisible nondivisible by r. We present two different constructions of an interior point of M, corresponding to the two cases.
Case 1: n (mod r)
We first construct a Up E Hrxr such that UoU 1, = 1,, and ||66 = K. Using the ansatz
Uo = (z u + y11" where x, y > O, we can find x and y by solving the system:
(UnUp 1, = 1,) 1 - (r 1)y² K/r (1(v = K)
x (r 1)y =1
The first equation gives x = (r 1)y. By substituting into the second equation, we obtain the following quadratic equation of y:
K = 0.
r(r 1)y² 2(r 1)y +1
By the quadratic formula and x. y > 0, we have the following solution
1 (1 ( 1)(K - 1)),
z 1)y
v 1)(K - - - -K-1
r(r-1)
Note that if K, still solve the quadratic equation, but without all positive solution. Denote = n/r, then U L, @ Ug an interior point of M,.
Case 2: (mod r)
Let us denote = |//r] and p = n mod r. We consider the ansatz
B1.1
B1
U [n
0] (1, B) for some block matrix B
Rrr
B
B
where
B1, = (x- v1,1
B1 = 21,1, r-pp
B2 (w 2)T =1 r-
B2,1 = yl,-p1,
19

---

Published as a conference paper at ICLR 2026
Additional to the constraints || K, and UUT = I,, assume that U1, c.In, L
ce ce and C. 1 for some G. and Ce. which are sufficient for UUT 1n Then we can find v, 2, and w by solving the system:
x ( (r p)z py - 1)2,
(q+1)x+ q(r - 1)y + (p- 1)y qu + q(r - 1) + ps,
(17)
( (r - p 1))) + 1)= + = 1
(1+1) (p- 1)²² (r- p)= + q(r pp) w² + pr-p-1);2) K.
The four equations correspond the following constraints, respectively: U1, C, In. U 1n co 1 crCe = 1, and |0, = K. From the first two equations, we can express x and in terms of z and (note that qr + p):
1
1
w
y =
1+
w.
(18)
n
72
By substituting (8) to the third and fourth equations of 17 we are left with system of quadratic equations of two variables:
01 + a2zw + azu² + c1 0,
(19)
b12² b2ow bzw? + 02 = 0,
where
nr qr pp-pppp-- + 1)/n,
02 = 2p(1 + 2q nn/n,
a3 -p(2q + 11/n.
bj (r + +p/n 1)(n 9).
b2 (1 pp/n -) + (( + -1).
by = g(1 p/n),
C1 K.
C2
Now our goal is to solve 19 By multiplying the first equation with bj and the second one with aj and subtraction, we can express 2 in terms of wt
aw + with = asdı anbs
b c 01
(20)
a162 a2b1"
a b2 a2bı
Suppose that w * O, we substitute (20) into the second equation and multiply by w The result is quartic equation:
(bia² + bzu + ba) + (2abbi + bb2 c2)w² 2 + 6bb, =
(21)
By solving
21 with aw + b/w 0, obtain
(2abb1 + 6b2 02) + /2abb) + bb2 02)² 4(a"b) ab2 + 63)6bbi
w
(22)
2(a2b; + ab2 + ba)
The proof is completed by combining (18) 20 and (22)
Proof of Theorem For "2," Uf [i, VJQ. Q. then UUT = 1,, + VVT and hence
UU 1 1 and (UU = K respectively, because 1, = and (r(V = K - For "C let U = PIQ denote the singular value decomposition with PT P = QQ = =, Since UU 1, = 1n,
the decomposition can be chosen so that PΣe; = 11/ = 1,, Sf PE ep), then V = (e,...., e,] = and | | - |E,, 1|2 --
In the final part, we first construct the inner approximation S of the tangent space
T(v.)) M ((V 1 (v, 0, QQT + +QQ 0} - = o. QQ [ h T
S.
h
0
We observe that
o
- h
i., v [i,
vl
Q = [V>
i
h
0
20

---

Published as a conference paper at ICLR 2026
and therefore the Jacobian operator is injective for all (V.Q)= E S:
p ØII" = || 1i vj v/ol
= [[lo vIl + [[v
in i || + |||l \22
/2 0v
where we used the fact that V = 0. Hence, the Jacobian operator is surjective, as claimed: dim (imainage (DV, (V, 2 dim(S) n(r 1) - 1 = dim(Tu MM).
E ADDITIONAL DETAILS ON RIEMANNIAN FORMULATION
E. EQUTVALENCE BETWEEN RIEMANNIAN OPTIMIZATION AND SQP
Let us consider constraint optimization problem.
min f(u),
"ER"
(23)
s.t. g:(u) = Vi.
where LICQ:
Σ (iV =0 - y: = AA.
holds for every in the constraint set. Then M M= E R": (.) - 0 Vi} is smooth embedded
manifold (Boumal 2023), and we may employ Riemmanian methods to solve 23
The Riemmanian method solves mingeM (u) by solving
min f(u) + (grad fu)), ú) + (Hess f(u) ú) WE M
||il|
(24)
each iteration. One can show that this is equivalent the SQP method that solve 23 by minimizing the Lagrangian
min f(), + (V.l(u, v ú) 1/2 V C(u, (")) ú) Lill
(25)
s.t. (Vgi(u), i)) Vi.
where
= ( + Σ vi g((u) and y (M) = arg minl L., v)]l-
Hence, our contribution is to efficiently solve the SQP subproblem by exploiting a block Alook--diagona. plus-low-rank structure in the Hessian, and the fact that there only r + rr+++++1 1)/2 < constraints. We provide more details in Appendix F.2 and Appendix E.3
To establish the equivalence, first observe that the search space of 24 and 25 the same. Indeed, the tangent space is T,, M {u: (V).). ú) = O Vi). Next, we write the expressions for the Riemannian gradient and Hessian (Boumal 2023 Prop. 3.61 and Cor. 5.16):
grad fu) Proj., (Vf(u)), Hess f(u)[t]] = Proj., (D, grad f(u)[tt),
(26)
where Proj.. (v) arg AT. M ||l ä||, and Dr. denotes the usual differential operator. We then obtain
grad (u) = Proj.. (Vʃ(u)) Vf(u) Σ M Vgl()) V ".(u,,y())
(27)
For the second equality, see Equation (7.75) in Boumal (2033
For the second order terms, we have (Hess [u] ú) = (V? [u)) ú) for all in T,, M by the facts that
21

---

Published as a conference paper at ICLR 2026
(i) Hess fu)][t]] = Proju (Vu
iii) PProj (v), (v, ui).
We obtain fact (i) by 266 and 27
Hess [u]][i] = Proj.
D grad [ [u)) (i)
= Projj, V2 (2 f(u)((())) [i]] Σ (Du v [gl(u) + Σ 0d v2 v2
EE
=Proj (V2 f(u) + Σ yi (u)
= -Proj Ww [(u, [u]).
Inside the projection operator, the term E, (Dy Vg vanishes because Ty M is a linear subspace, and (T,, M)) span (Vgi Fact (ii) is due that the projection operator is self- adjoint and that the projection of any tangent vector is itself. For more on the connection between SQP and Riemannian Newton method, we refer to Absil et al. (2009); Mishra & Sepulchre (2016).
E.2 GENERAL FORM THE RIEMANNIAN GRADIENT AND HESSIAN
Throughout this section, we a bar over a function defined on the manifold M to denote its smooth extension defined on a neighborhood of M so that the Euclidean gradient and Hessian can be defined on MM. Namely, for MM R" use f. (NM) - RT to denote the smooth extension. The notations grad f and Hess f denote the Riemmanian gradient and Hessian; V and 2 denote the Euclidean gradient and Hessian.
For manifolds that can be defined by
min f(U),
M = (U E Rax A(UUT) + B(U) c),
DEM
Rm linear operators, and R its tangent space can where A: ROxx R and B: Raxr be written as:
To (Ú _axr AQÜUT +UÚT)+ B(Ü) 0}. We call the function A(UUT) + B(U)= = as the defining function of M. Let us denote
L(Ú) == A@UUT )
(28)
Immediately, we see that Ty M = ker L). and the adjoint operator
L (y(r) 2 A (y)U B (y). The projection operator onto the tangent space is defined to be
(29)
Proju (W) = arg min |W Ü|| ÚETo M
We also know that ker(Z) = imago(L") L*). Thus, by the orthogonal projection theorem, we may that
Proje (W) W w W - L
(30)
where arg min,, IW-LL" (y) is the solution the linear system W L')) E ker(L) Tu M. and both Proju (W) and ÿ are unique.
Consequently, the Riemannian gradient and the Hessian matrix-vector product have the following form:
grad f(U) = Projo (VJ(U)) Vi(U) +(A (vo )U + +B (yc)) Hess f(U)[U] = Proju (D grad ((U)((]) Proi. (v2 [(U))( 2[A (vo))(),
(31)
where yu -Ø in (30) the unique Lagrange multipliers
Vf(U) + 2 [ + B (T)
y = arg min m
(32)
For detailed proof, see Boumal al. 2020
22

---

Published as a conference paper at ICLR 2026
E.3 EFFICIENT COMPUTATION OF THE RIEMANNIAN GRADIENT AND HESSIAN
We first write down the Euclidean gradient and Hessian for our objective function, and then explain how compute the Riemannian counterparts efficiently. Specifically, we show that yo can be computed in O(r + r²) time.
We decompose the objective function as f f f1 + uff, where
log Pi, (V, fiV. = (C,VVT), Q2 Σ
3.j
For fI, the Euclidean gradients are Vvfi(V.Q) = 2CV and VQfi(V,Q) = -0, respectively. For f2. the gradients are
= -V
where U== == QV, and U(-1) the element-wise reciprocal, i.e., [U(-1)) = UTJ- For conve- nience, also define Q DQ with D 10,-1 I--1] Note that Q is simply Q without the first row. Putting these together, we obtain
Gv = Vv.( = 2CV - pU (T)
G =-N
The Euclidean Hessians can be given in vectorized form as
21 1 60 C
0
dvec(U [Jv
Jql
0
Hvq
Ja
Hov
0
where
Jv = Q @ In
JQ = 1, 2 QU
and
HvQ = OD K(m)) = =HQv
with U(-2) being the element- square of U(-1), dvec(.) diag ddag[vec(-)) (-1 and (( ((())
denoting the commutation matrices (Magnus & Neudecker 2019, Sec. 3.7). We can then compute
the Riemannian gradient and Hessian-vector-poodne according to 31 and (2). In the remainder of this section, show how to efficiently solve (32
The manifold consider, M == Orth(r), can be written
M = { (V.Q) (r R 1 V 0, 0, K svec 1,) - =0] =
where svee denotes the symmetric vectorization (Klerk, 2006 Appx. E). Note that there are no cross terms (VQT or QVT) in the defining functions of M.. Thus, we can treat the defining functions with respect V and Q separately. The corresponding terms are
Av(VVT)= = [ 0 In(Vitt]
BV( = 1 0
Aq(QQT) = svec(QQT)
and
0
cV ==
K
and cQ = svec(I,).
Mimicking
(88 and 29
we use the notation Lv, Lo Ly, and L2 respectively. For any (V1,
Rr-11 and ys E 11/22 we have
L (Vi y2) 1, 2 V and L2(ys) smat(ys)
where smat is the inverse of svec, that is, smat(svec(MM) (M)) = M for all symmetric matrices M. By solving the linear systems
Gy (v1, ker(Tv)) and Go I. (ya) ker(T.q))
23

---

Published as a conference paper at ICLR 2026
obtain the following closed form solutions:
51 1 G Im j/2
V) 2(K 1)
1 5a svec(GoQ QQ))
(33)
The computation of D1, 82, and y in total requires (nr + + time. Therefore, can compute gU - (3), 22, V3) with the same cost. Given a Euclidean gradient and a Euclidean Hessian-vector product, we may write out explicitly the Riemannian gradient:
[ Tm
1 1
(Gv K-1
GQQT-QGy -QG
grad f(V.Q) =
Gv
(34)
2
and the Riemannian Hessian-vector product:
Projv (2S(V, 21v (Gv, VGr.VV
K-1 GoQ + +QCae
Hess S(V,
Projo V"
2
E.4 FEASIBLE INITIAL POINT
K is necessary and sufficient for the existence of an interior point of M. The following Lemma shows the necessity of > K. When K, the structure of the unique U Ra*k in Lemma can be explicitly written as
In this section, we show that
VIG G.
1a
1xx
(35)
VIGal
VIGK
where 1G, E {0,11 1}"h denotes the binary vector with its support being Gk.
Lemma 5 (Isolated feasibility when r = K). Let -M RT and M4 = M'
where RT*K (U E Ba*K U 0}. Then, we have: (i) M+ M4; (ū) M+ and only U J is a group assignment matrix defined in (35): (iii if is a group assignment matrix, then the intersection of the tangent space Ty M and the cone Ct = {V E Raak :. 0, Vu.g 1 0} is trivial, i.e.,
To MrMMnC nC: {0).
In Lemma M when r > K. Here, we present the construction when n = qr is an integer multiple of r. Let Un = (x y) yy,1 where x = (1 v(r 1)(K 1)) and y = r-1(vr - / K 1).
moreover provide complete analytical construction for an interior point of
Then U i, @ Uo is an interior point of M. Next, we show how to compute the pair V, Q corresponding to the interior point U by SVD. For given U E Mr, let U = Pu be the SVD of U. We find (V,Q) such that U VQ by V sgn(Puran))PPo))... and Q aT
E.5 LIPSCHITZ CONTINUITY OF PENALTY
To apply the guarantees in Section 2.3
need to take of the logarithmic penalty in (3) since it does not have Lipschitz gradients nor Hessians over its whole domain. The standard workaround, widely used in the analysis of nonlinear interior-point methods, is observe that all iterates U₂ (V&,, Q&) remain strictly feasible. Consequently, the penalty could be modified by Huber-style smoothing, where d min,,,,((UR)., > O:
log
2 8
r(x)
log 8 + () (r (x-
(x
5)3
28
28
The function r(x) both concave and has Lipschitz Hessians. Therefore, the guarantees in Section apply. The smoothing is only needed for theoretical purposes. In practice, apply the Riemannian algorithms directly to log z, and not to ræx). Since we have assumed that all queries satisfy x 2 the actual behavior remains consistent with the smoothed model.
24

---

Published as a conference paper at ICLR 2026
F EFFICIENT IMPLEMENTATION AND COST OF BISECTION SEARCH
To implement the proposed method, we vectorize the input as = (e, q), with v = vec (oo(VT) and
q vec(Q). Since -XX we rewrite the cost function form Appendix E.3 as
f(u) = -1| || VII2 n1, log (VVV
and define the constraint functions
(ı(u) = 1 V.
g2(u) = | - (K 1),
ga(u) = svec(QQ I.).
The Jacobian J and Hessian H are computed analytically, as Appendix E.3 their sparsity.
in order to exploit
Since have restricted ourselves to optimizing over p T,, M,, may replace the Riemannian Hessian H with H Cuu, ya). The key observation is that H admits block-diagonal-puss low-rank structure.
For convenience, we list some of the derivatives in this section. The (Euclidean) Jacobian of the constraints can be written in block form as
J 1 0
=
J 2n
0
J3q_
where
J10 = (Tr-1
11,1 (4-1
J2v 2v,
J3g S,(1 + (Q @ I,),
where S, E (R((++1)/2x" the matrix satisfying S vec(Q) = SQ)
Computing the second-order derivatives of @1: 92 is straightforward, since J1o is constant in V and Jz. is linear. To compute the second-ordes derivatives of gs. since we only need to compute Ej4r((41)2) Ma.j V (3.ss(()) we note that any y E Rr(r+11/2
Jzgy = vec(smat(y))))
Consequently, we have for any y Rr((+1)/22
Σ
vV 93.j 2I, @ sya
jelr(r+1)/2
Collecting results, we have
How BBB
Hog Hgg
Hgv
with
= V2(X 8 I--1)
(36)
Hos = nl, & Q) dvec U(-2) (Ta 8 (T) +j27,
(37)
Hgg (T, & V dvec U (2) (T. v) 2I, smat(ūa).
(38)
and
Hug -u(U 8 D) + u(1, @ Q) (mr) dvec(U
HTT
(39)
where V2, 33 follow (33
To solve the saddle point problem arising from subproblem linear system
14
via bisection search, we solve the
H XI BBT
Hog
J
Hqy J
Hay
Al
J
=
fg
Jn
r
25

---

Published as a conference paper at ICLR 2026
Repartitioning along the lines yields:
Kn
K2
b1 b2
K21
K22
32
with æ1, bj E [=(r-1) and 22, b2 E *m*+m In particular, we observe that the block K11 has the form K11 D11 BB", where D11 Hev + XI is block-diagonal, with n blocks of and B has at most dr columns. Therefore, we instead solve
D11
K2
b
K
K22
=
b2
B
I
First, it costs n(r 1) = O(nr*) time to invert D1I- Afterwards, forming and solving the size m + x² + rd Schur complement problem:
(L22 D1 L12)
q
b2 [
L D
bi.
(40)
2
where
K22 0 0
L12 == [K21
BJ
L2 =
I
cost (d 6 ) time. In the end, we substitute to recover
DTi b1
12
in (d++ time, and apply retractions to and y. In total, it takes ( + + 6 + rdd3) time solve the system, which is indeed n poly(r, d). Putting pieces together, pseudo-code of our Riemannian method is shown in Algorithm
26

---

Published as a conference paper at ICLR 2026
G ADDITIONAL NUMERICAL RESULTS
We collect additional numerical results this section.
Dataset visualization. Figure g and Figure 2 display the first two principal components of the GMM dataset and CyTOF dataset, respectively.
y= 2
y= 0.5
7.5
5.0
2.5
-2
-25
5.0
-7.5
-6
5
-5.0
-2.5
0.0
2.5
5.0
Figure 8: Visualizing the effect of the separation parameter in GMMs. As y decreases, the clusters become increasingly difficult to distinguish.
4
2
0
-2
-5.0
-2.5
0.0
2.5
Figure 9: Visualization of the CyTOF dataset. Two clusters exhibit significant overlap, implying the difficulty of clustering.
Benchmark on CIFARI0. Figure TO shows the performance of our method on the CIFARIO image classfication dataset. The original 32 32 color images were processed using pre-trained Vision Transformer (ViT-B-16 16) (Dosovitskiy et al. 202T) to extract 1i8-dimensional features, which were then reduced to d 50 using PCA. We setting. From the dataset, 25000 images across five classes were selected; In each trial, we draw subsample of 1,000 images and perform clustering. We repeated this procedure 50 times.
pre-trained model avoid breaching the unsupervised
Generalization to kernelized K-mean. Since the data enter the problem through the Gram matrix C = T in Equation 33. extending the implementation to the kernel K K-means is straightforward.
27

---

Published as a conference paper at ICLR 2026
Mis-clustering Error
0.4
10
8
0.3
5
0.2
0.1
Ours
NLR
SC
KM
10
0
5
Figure 10: Real-world benchmark on the CIFARIO data. (Left): comparison to other state-of-the-ane methods. (Right): first two principal components of the image embeddings.
We adapted our implementation to accept a kernel matrix as input and conducted experiments using the RBF kernel toy datasets. The results shown Figure
Estimated Clasters
Global Suboptimality
Grudient Norm
1.0
Sa #
10
10'
10-3
10-2
10-6
10-9
10
10-12
10
250
300
750
250
500
750
1.0
10"
100
10-3
0.5
10
10-6
10"
10-12
10 6
-055
Iterations
Figure 11: Clustering with the RBF kernel. Demonstrated on the Circles and Moon datasets from scikit-learn. The kernel bandwidths are and 15, respectively.
Nevertheless, without access to the factored matrix, described in Appendix F and therefore we lose the linear-time guarantee for solving the subproblem. To extend the linear-time core claim would further require low-rank approximations of kernel matrices C, e.g. via Nyström. We leave this as important future work outside the scope of this particular paper.
longer perform the sparse partition
Impact imbalanced clusters. Figure 12 shows the impact imbalanced cluster sizes. this example, the clusters contain approximately 10%, 20%, 20%, and 50% of the 1000 data points. We compared our method with NLR and spectral clustering (SC), omitting K K-means++ due to its substantially higher error, which would distort the scale of the results. As the imbalance among clusters increases, we noted that the convergence becomes less stable with the default parameters, so increase the rank from K to K + 3. The new imbalanced cluster size experiment further
28

---

Published as a conference paper at ICLR 2026
reveals two things: (i) the performance of all methods under comparison degraded, and (ii) method is the most robust (with the narrowest error distribution) among all.
Mis-clustering Error
Global Suboptimality
Gradient Norm
10²
10a
0.006
10'
10
0.004
10
10-5
10
0.002
10-9
10
0.000
Ours
NE.R
SC
100
200
300
100
200
300
Iterations
Itcrations
Figure 12: Clustering performance and convergence behavior when cluster sizes are imbalanced.
To explain this empirical observation, it is known that similar (and more general) recovery threshold in the motivating K-means SDP to Equation (6) holds for unbalanced cluster sizes. In particular, the exact recovery threshold depends the minimum of harmonic means of the cluster sizes, ie. m = minisk mail:sss//((2mmm (2ngn))/(ng - ni), in the following way (Chen & Yang 2021b):
d
θ 2 1+
m log n
In the special case K equal clusters, m n/K attains the maximum possible value. Hence, when the cluster imbalance is significant, a larger centroid separation is required for the proposed algorithm to achieve clustering accuracy and convergence behavior comparable that observed in the equal-cluster case.
Runtime scaling with respect to r and Figure shows the average per-iteration runtime as as r and vary. The setting is the same that of Figure 2 with 1000 fixed. The observed scalings for r and d are roughly (rr*) and (dd), consistent with the leading terms of ((r + +r*d Because dominates leading constant, the contributions of the 5 and d terms are unlikely become significant until r and d are extremely large.
Elapsed Time (second)
0.15
0.015
0.10
0.010
0.05
0.005
0.00
6
8
10
12
20
40
60
80
100
d
Figure 13: Average per-iteration runtime versus rank and dimension, 1000. The measured runtimes exhibit approximately (rr3) and (dd) scaling (dashed lines).
Robustness to initialization. As illustrated by Figure I our method is robust to initialization, all 50 trials successfully converged to second-order optimal solutions. Although the solutions differ
29

---

Published as a conference paper at ICLR 2026
(Figure 14), their corresponding membership matrices Z are close to each other (Figure 15), and yield identical clustering result. Moreover, the minimum eigenvalues upon convergence form distinct clusters that align with clusters in the recovered membership matrix Z, shown in Figure [5] These local critical points consistently produce perfect clustering, indicating that they remain close to the global optimum.
0.10
0.08
18
0.00
Figure 14: Difference between the solutions. First 25 rows of selected solution U obtained from the
global optimality experiment described in Section 4
Frohenius Distance
6 3
0.016
18 15 12 9
0.014
0.012
2
24
30 27
0.008
48 45 42 39 36 33
0.006
0.004
0.002
0.000
15 21
27 30
36
42 48
Figure 15: Similarities of the membership matrices. Frobenius distances between the membership matrices obtained from the global optimality experiment in Section sorted according to their corresponding minimum Hessian eigenvalues.
Comparison with another Riemannian elustering method. Inexact Accelerated Manifold Proxi- mal Gradient Method (I- (I-AManPG) by Huang et al. (2025) is a recent first-order Riemannian method for solving general problems of the form min
f(X) +X|X),
(41)
s.t. X F = =x x span(X))
We evaluated its performance on the clustering problem using the on the CyTOF dataset with 50 repetitions (same setting in Figure 3. The results shown in Figure 16 While I-AManPG is
30

---

Published as a conference paper at ICLR 2026
generally fast and accurate, its median error is higher than that of our methods. In particular, several runs of I-AManPG exhibited large errors, indicating convergence failures.
Mis-clustering Error
0.4
0.3
0.2
0000
0.1
Ours
I-AManPG
Figure 16: Comparison with I-AManPG using CyTOF. Performance of I-AManPG is comparable to other clustering methods. However, it suffers from convergence failures from time to time and requires careful tuning. Our method again demonstrated its accuracy and stability.
Hyperparameters tuning. the various numerical experiments, we observed that a smaller value of y led to more accurate solutions but at the cost of slower convergence. Therefore, recommend selecting the largest possible y that does not trigger the phase transition. A good heuristic we found to choose such that the initial penalty term μf2 remains less than 20 times the main term fi in the loss function. The onset of phase transition is also easy to notice, as the algorithm will quickly stagnate and terminate in just few iterations. If higher accuracy is desired, one reduce ja gradually, using the solution obtained with a larger as initialization. This warm-start strategy significantly speeds up convergence compared to using small from the start.
As noted earlier, increasing can also improve accuracy, likely because it improves the problem landscape. However, due to the poly(r, (r, d) runtime scaling, caution must be taken when deciding whether to decrease strategy: begin by decreasing l with K 1; if this does not produce satisfactory convergence or leads to excessive computation time, try increasing r to K + 2.... K + for some small c, and repeat the y-tuning process.
increase If tuning (μ, together is desired, we suggest the following
The other hyperparameters in A D influence the speed of the inner optimization. The initial multiplier À affects only the number of inner steps required during the first iteration. We recommend doing simple trial run with only two iterations; the resulting optimization history typically offers a reliable guide for choosing an appropriate initial scale for For the other two parameters, suggest setting slightly smaller than K4. Empirically, we found x_ 1.1 and 64 = 1.3 work well.
Additional convergence plots. Figure 17 Figure 18 and Figure 19 illustraat the convergence of
method GMM with different parameters and on the CyTOF dataset, demonstrating its stability across different datasets.
Hardware information. All experiments in this work were conducted on a machine equipped with single Intel Core 1490KK CPU and 32 GB of RAM.
31

---

Published as a conference paper at ICLR 2026
Global Suboptimality
Gradient Norm
Minimum Hessian Eigenvalue
109
2000
10
10a
10
10a
-2000
10-1
10
-4000
6000
Iterations
Iterations
Iterations
Figure 17: Convergence our method on GMM data with perfect separation (n = 500, = 1.0). The loss value steadily decreases iterations and converges rapidly near the optimal point. This example achieved a perfect final clustering result the end.
Global Suboptimality
Grudient Norm
Minimum Hessian Eigenvalue
10a
10
10
10
10-1
10"
-1000
10-1
10-8
-1500
2000
0
100
200
200
100
200
lterations
Iterations
Heratiioss
Figure 18: Convergence of our method on synthetic Gaussian mixture data with low separation (n 500, 7 0.25).
Global Suboptimality
Gradient Norm
Minimum Hessian Eigenvalue
10
10'
10a
10-1
10
-000
10-
101
-4000
10-9
10-3
6000
100
100
0
100
200
Iterations
Iterations
Iterations
Figure 19: Convergence behavior of our method on the CyTOF dataset.
32

---

Published as a conference paper at ICLR 2026
License of assets. The MANOPT solver is distributed under the terms of the GPLv3 license: the PYMANOPT solver is released under the 3-Clause BSD license; the CyTOF dataset is the work of Levine et al. 2015), cleaned and distributed by Weber (2015) under the MIT license.
Code. A Python demonstration of the proposed method is available at https:///itte com/ Francis- fisu/kmean manifold This repository provides example scripts for the GMM and
CyTOF experiments presented in this work.
H
PSEUDOCODE
This section presents the pseudocode of the proposed method. For its derivation, see Appendix F
Algorithm
Riemannian Second-order Method
Require:
Data X,
Initial point (Vo, Qo).
Initial multiplier Increment/decrement factors (K4,K-),
log-barrier penalty μ.,
Max number of outer/inner iterations T and B.
1: (V, < (V, Qo) for i. do
2:
3:
Vectorize the input: u + |vec(V vec(Q)]T;
4:
Compute the current loss: L f(u);
5:
Compute the Riemannian gradient G + Vf(uu;
6:
Compute the Jacobian J < (());;
7:
Compute H V C(u, Vu) as in Appendix F
8:
for 1,......... ..... B do
9:
0, < SOLVEINNER G,
10:
Reconstruct V, Q from vector ú g
11:
(V.Q) + (V, Q)
12:
Compute the new loss L' from (V, Q)
if L' then
14:
À < A/K-
break
16:
else
17:
4 A K.
18:
output (V,Q)
19:
20:
function NNR, J.
Add A/ the V and Q blocks of H:
22:
Form block matrices as in (40):
23:
S
L22
Di L12
b2
L
D bj
0
D
b
L
12
7
26:
return
33


Reviews:
[{"review_id": "aRBG10QkFx", "reviewer": "Reviewer_9Enk", "summary": "This paper proposes a scalable second-order Riemannian optimization algorithm for solving the K-means clustering problem through a low-rank semidefinite relaxation. The authors reformulate the nonconvex problem as a smooth unconstrained optimization problem over a product manifold and apply a cubic-regularized Riemannian Newton method. The algorithm is claimed to achieve linear per-iteration complexity with provable convergence to second-order critical points.  \n\nHowever, several issues remain. The key assumption (Assumption 1) seems inconsistent and lacks justification; the algorithmic novelty is limited.   If the authors can address my concern, I would raise my score.", "questions": "Why is the cubic-regularized Newton framework preferred? The original method (Agarwal et al., 2021) solves subproblems via *Lanczos iterations*, whereas this paper adopts bisection search. The authors should discuss those two approaches.", "limitations": "", "rating": 4, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "- The mapping from the constraint manifold $\\mathcal{M}$ to the product manifold provides a principled and novel way to represent K-means as a smooth manifold problem.   This reformulation reduces the projection cost from $O(n^2)$ to $O(nr + r^3)$.  \n- By exploiting block-diagonal-plus-low-rank Hessian structure (Appendix E), each iteration scales linearly in \\(n\\) in the cubic-regularized Newton approach.  \n- Use bisection search to solve the subproblem in the cubic-regularized Newton approach.", "weaknesses": "- Assumption 1 is difficult to interpret. The problem (1)  does not require \\(U > 0\\),  (U > 0\\) enforced in problem (3), which is only a sufficient but not necessary condition for $Z = UU^\\top > 0$. Hence, the assumption is logically inconsistent with the model definition.  There exists the case that the optimal solution $U_{ij}<0$ .  In addition, this assumption is verified only empirically (Fig. 1), with no analytical justification or example provided.  Since all global optimality claims rely on this assumption, its ambiguity weakens the theoretical contribution.\n\n- The method essentially applies a standard cubic-regularized Newton algorithm to a reformulated manifold problem, which not  introduce new algorithm, or adaptive strategies, or theoretical improvements.\n\n- Although the reformulation simplifies the constraint structure, computation of gradients and Hessians remains costly.  I recommend the authors discuss hybrid methods that separate simple manifold constraints (handled via Riemannian optimization) and other remained constraints (handled via augmented Lagrangian function).  Related work worth discussing includes:  \n\n  - Wang, J., & Hu, L. (2025). Solving low-rank semidefinite programs via manifold optimization. Journal of Scientific Computing, 104(1), 33. \n\n  - Monteiro, R. D., Sujanani, A., & Cifuentes, D. (2024). A low-rank augmented Lagrangian method for large-scale semidefinite programming based on a hybrid convex-nonconvex approach. arXiv preprint arXiv:2401.12490. \n  -  Hou, D., Tang, T., & Toh, K. C. (2025). A low-rank augmented Lagrangian method for doubly nonnegative relaxations of mixed-binary quadratic programs. Operations Research.  \n  -  Wang, Y., Deng, K., Liu, H., & Wen, Z. (2023). A decomposition augmented lagrangian method for low-rank semidefinite programming. SIAM Journal on Optimization, 33(3), 1361-1390.\n\n- Important components—algorithm pseudocode, subproblem solvers, and retraction details—are placed only in the appendix.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "8cj9dyrGRF", "reviewer": "Reviewer_9AoF", "summary": "This paper tackles the nonconvex low-rank factorization of SDP relaxation of K-means clustering. By penalizing the nonnegative constraints, the authors view the resulting problem as a manifold optimization problem and apply the Riemannian cubic regularized Newton to obtain a 2nd-order solution.", "questions": "1. Is the part of lines 311-355 new result? If this is not new, it would be better to cite some references.\n2. Could the authors explicitly explain why Riemannian cubic regularized Newton can overcome the ill-conditioning and RGD/RTR cannot?", "limitations": "", "rating": 6, "confidence": 5, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "1. Develop a fast algorithm from manifold optimization viewpoint (Theorem 2), which is novel at least to me.\n2. The subproblems of Riemannian cubic regularized Newton can be solved efficiently.\n3. Empirically, 2nd-order critical points are globally optimal.", "weaknesses": "1. This paper only provides sublinear convergence rate. Is it possible to prove stronger local convergence rate results? For example, quadratic/superlinear convergence rate. Fig1 shows fast convergence rate in practice.\n2. The benign nonconvexity is described in Assumption 1. Is it possible to say something about benign landscape? For example, related work of benign nonconvexity results.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "6NsW6V1QZo", "reviewer": "Reviewer_dh1D", "summary": "The paper proposes a new smooth, unconstrained Riemannian reformulation of a nonnegative low-rank factorization of the K-means SDP, enabling second-order optimization with feasibility preserved via retractions. The key idea is to map the original constraint set $\\mathcal{M}$ onto a product manifold $\\tilde{\\mathcal{M}} = \\mathcal{V}\\times \\mathrm{Orth}(r)$ via a submersion, which yields simple $O(nr + r^3)$ retractions and efficient expressions for Riemannian gradients/Hessians.\n\nBuilding on this geometry, the authors design a cubic-regularized Riemannian Newton method that solves each Newton subproblem in $n\\cdot\\mathrm{poly}(r,d)$ time by exploiting a block-diagonal-plus-low-rank structure and a bisection scheme on the regularization parameter. The overall algorithm finds $\\epsilon$–second-order points in $O(n\\cdot\\epsilon^{−3/2} \\cdot\\mathrm{poly}(r, d))$ time.\n\nThe work hinges on an empirical \"benign nonconvexity\" assumption: in regimes where the convex $K$-means SDP recovers ground truth, all approximate second-order critical points of the nonnegative low-rank model are near-global optima. The authors provide extensive empirical evidence for this behavior.\n\nExperiments on synthetic GMMs and CyTOF data show that the proposed second-order method converges in hundreds of iterations, achieves similar or better clustering accuracy than the strongest first-order baseline (NLR), and substantially reduces time despite costlier iterations. It also outperforms prior Riemannian $K$-means methods and classical RTR/CG solvers that struggle with the log-barrier’s ill-conditioning.\n\nThe paper supplies detailed derivations: manifold geometry, LICQ verification, submersion proof, closed-form Lagrange multipliers for projections, efficient Hessian-vector products, feasible initialization (and necessity of $r > K$), and implementation details for the linear-time inner solves.", "questions": "1. Can you provide partial theory toward Assumption 1? For example, under the separation in Chen & Yang (2021b), and mild overparameterization $r = K + O(1)$, can you show absence of spurious second-order points in a neighborhood of the ground-truth factor, or establish a strict saddle property?\n\n2. How robust is the method to mis-specified $K$ and to cluster imbalance? Could you include experiments where $K$ is under/over-estimated, and where cluster sizes vary significantly, and report both accuracy and convergence behavior?\n\n3. Could you explore alternative penalties that retain positivity while improving conditioning (e.g., softplus smoothing, additive offsets, or barrier homotopy/continuation schedules), and compare convergence to your log-barrier?\n\n4. What is the practical guidance for selecting $r$ beyond $r = K + 1$? Are there cases where larger $r$ improves optimization (escaping poor basins) or statistical robustness, and what is the runtime trade-off empirically?\n\n5. Can you extend experiments to kernel $K$-means and non-Gaussian mixtures, or real vision/NLP datasets where the Gram matrix is built from learned features? This would help demonstrate the generality implied in Appendix A and the stated manifold framework.", "limitations": "", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "Conceptual advance: The submersion to a product manifold with simple retractions is elegant and practically impactful, removing the main bottleneck that hindered prior Riemannian approaches to $K$-means (expensive retractions and feasibility maintenance).\n\nAlgorithmic engineering: The cubic-regularized Riemannian Newton solver is carefully tailored—analytical gradients/Hessians, efficient tangent projections, and a block-diagonal-plus-low-rank exploitation for the inner linear systems. The bisection-based solver for the regularization parameter is simple and reliable.\n\nClear bridge to theory: The paper situates the contribution within the exact-recovery phase transition for the $K$-means SDP and explains how second-order points suffice under the benign nonconvexity hypothesis. The LICQ and manifold calculus are handled rigorously; the smoothing argument for the log penalty clarifies the use of second-order guarantees.\n\nStrong empirical evidence: On both synthetic and real data, the method exhibits rapid, stable convergence to second-order points with near-zero mis-clustering; it consistently reduces iteration counts by orders of magnitude compared to NLR and RTR, translating to 2-4x faster runtimes despite more expensive steps.\n\nReproducibility and completeness: The paper provides explicit formulas, complexity accounting, initialization constructions, hyperparameter guidance, and ablations (rank $r$, penalty $\\mu$, comparisons with several baselines), which make the contribution actionable for practitioners.", "weaknesses": "1. Assumption 1 (benign nonconvexity) is purely empirical in this work. While the authors motivate it with analogies to Burer–Monteiro and back it with experiments, the lack of any partial theoretical characterization (e.g., under separation/noise conditions and mild overparameterization $r > K$) limits the scope of the main claim in average-case regimes.\n\n2. Dependence on the log-barrier: Although handled well algorithmically, the severe ill-conditioning induced by the log term drives both the design choices and some limitations (e.g., RTR/CG underperform). A discussion or experiment on alternative barrier/penalty designs (e.g., smooth-plus hard positivity via projections on $U$) could strengthen the case or show robustness.\n\n3. Sensitivity to $\\mu$ and feasibility interior: The method requires strictly positive $U$ and shows a phase transition when $\\mu$ is too large. While the paper provides heuristics, a more systematic procedure (or adaptive schedule with safeguards) would make the solver more turnkey across datasets; also, the necessity of $r > K$ to ensure interior feasibility may be restrictive in memory-limited settings.\n\n4. Scalability constants: The per-iteration complexity is linear in $n$ with $\\mathrm{poly}(r, d)$ factors; however, the inner solves involve $r^3$ and $d$-dependent Schur complements. Scaling with varying $r$ and $d$ would help practitioners understand the limits and guide parameter choices.\n\n5. Generality beyond GMM: Although the manifold formulation applies to kernelized $K$-means, empirical validation is focused on GMMs and one CyTOF setup. Broader tests (imbalanced/many clusters, higher $d$, other real datasets, kernels) would better support claims of robustness.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "YPt43MgKce", "reviewer": "Reviewer_dYYy", "summary": "The paper tackles the computational difficulty of solving the nonconvex, constrained formulation of K-means clustering. It reformulates the low-rank SDP relaxation as a smooth optimization over a Riemannian manifold, so that constraint feasibility is handled intrinsically. This allows the use of a cubic-regularized second-order Riemannian Newton method that guarantees convergence to second-order critical points—empirically aligned with global optima in K-means.\n\nTo make the approach scalable, the authors introduce a product-manifold factorization and exploit the Hessian’s structure to achieve linear-time per-iteration complexity. Experiments on Gaussian mixture and real-world datasets show faster convergence and lower mis-clustering rates than state-of-the-art first-order methods, demonstrating both theoretical soundness and strong empirical performance.", "questions": "1. In line 295, how sensitive is the algorithm to the initial point $(V_0,Q_0)$?\n\n2. Could this approach generalize to kernelized or probabilistic K-means variants? If so, how would the manifold and retraction structures adapt?\n\n3. How is the second-order condition in Eq.(10) verified numerically? Does this influence empirical runtime?\n\n4. As discussed in Appendix C.1, the Riemannian cubic-regularized Newton is equivalent to certain SQP methods. What are the computational advantages compared to projected Newton or augmented Lagrangian solvers?\n\n5. In the proof of Lemma 4, the fourth equality is incorrect, but I believe it's just a typo and not a significant issue.\n\n6. In Eq. (11), is the efficiency of the Hessian 1/2?\n\n7. Typo on line 1127: this is not a quintic equation, but a quartic one.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "1. The paper successfully integrates second-order Riemannian methods, specifically the cubic-regularized Newton approach, into the K-means SDP framework. Moreover, The proposed product-manifold factorization resolves the computational bottleneck of expensive retraction operators in previous works.\n\n2. Demonstrates consistent improvements in convergence speed and clustering accuracy over leading baselines, including NLR, RTR.\n\n3. Sections 2 and 3 provides a concise review of SDP relaxations and information-theoretic thresholds. Theorem 1 and Theorem 2 are well-chosen to connect existing Riemannian convergence theory to this specific manifold formulation.\n\n4. Reduces per-iteration complexity to linear in n, making second-order methods computationally competitive. The approach can be generalized to other manifold-constrained ML problems.", "weaknesses": "1. Assumption 1 (Benign Nonconvexity) remains unproven. The paper heavily relies on this assumption for its theoretical motivation but provides only empirical justification. A rigorous explanations or partial theoretical evidence would substantially strengthen the work.\n\n2. In the numerical experiments, only a single real-world dataset (CyTOF) is used. Additionally, there is no ablation study on manifold dimension, sensitivity to initialization, or robustness to cluster imbalance.\n\n3. The main theorems of complexity are adaptations of known Riemannian results rather than new convergence analyses specific to K-means.\n\n4. The paper’s contribution seems somewhat incremental. It reformulates the K-means problem and employs a standard cubic-regularized Newton method, with the primary advancement lying in a more efficient computation of the Newton step by exploiting its structure.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)