Riemannian Optimization on Relaxed Indicator Matrix Manifold
Authors: Jh Yuan
Fangyuan Xie
Feiping Nie
Xuelong Li
Venue: ICLR 2026 Poster (2026)
Keywords: Optimization
Clustering
Graph Cut
Primary area: unsupervised, self-supervised, semi-supervised, and supervised representation learning

Abstract: The indicator matrix plays an important role in machine learning, but optimizing it is an NP-hard problem. We propose a new relaxation of the indicator matrix and compared with other existing relaxations, it can flexibly incorporate class information. We prove that this relaxation forms a manifold, which we call the Relaxed Indicator Matrix Manifold (RIM manifold). Based on Riemannian geometry, we develop a Riemannian toolbox for optimization on the RIM manifold. Specifically, we provide several methods of Retraction, including a fast Retraction method to obtain geodesics. We point out that the RIM manifold is a generalization of the double stochastic manifold, and it is much faster than existing methods on the double stochastic manifold, which has a complexity of \( \mathcal{O}(n^3) \), while RIM manifold optimization is \( \mathcal{O}(n) \) and often yields better results. We conducted extensive experiments, including image denoising, with millions of variables to support our conclusion, and applied the RIM manifold to Ratio Cut, we provide a rigorous convergence proof and achieve clustering results that outperform the state-of-the-art methods. Our Code is presented in Appendix H.

Published as a conference paper at ICLR 2026
RIEMANNIAN OPTIMIZATION ON RELAXED INDICATOR MATRIX MANIFOLD
Jinghui Yuan, Fangyuan Xie, Feiping Nie Xuelong Li School of Artificial Intelligence, Optics and Electronics (iOPEN), Northwestern Polytechnical University, Xi'an 710072, P.R. China
yuan jhéma 11. nwpu edu. - en, xi xie faggaanyanemaaa .ww.. en deipingnie@gmall com; 1 li wwwppuuduu.... on
ABSTRACT
The indicator matrix plays an important role in machine learning, but optimizing it is an NP-hard problem. We propose new relaxation the indicator matrix and compared with other existing relaxations, it can flexibly incorporate class informa- tion. We prove that this relaxation forms a manifold, which we call the Relaxed Indicator Matrix Manifold (RIM manifold). Based on Riemannian geometry, we develop Riemannian toolbox for optimization on the RIM manifold. Specifically, we provide several methods of Retraction, including a fast Retraction method to obtain geodesics. We point out that the RIM manifold is a generalization of the double stochastic manifold, and it is much faster than existing methods on the double stochastic manifold, which has a complexity of (n²)), while RIM mani- fold optimization is (()) and often yields better results. We conducted extensive experiments, including image denoising. with millions of variables to support our conclusion, and applied the RIM manifold to Ratio Cut, we provide a rigorous convergence proof and achieve clustering results that outperform the state-of sate-of-thhe-ar the- art methods. Our Code is presented Appendix H.
INTRODUCTION
Indicator matrices play crucial role in machine learning (Mo al., 2025: Li et al., 2024a; Tsit- sulin et al., 2023), particularly in tasks such as clustering (Fan et al., 2022; Macgregor, 2024) and classification (Shi et al., 2024). For a problem with samples and e classes, the indicator matrix F E Ind where Ind = {X E "R Xi E {0, X1, - 1 and 1. is the column vector of ones of size c. The optimization of indicator matrices, which can be seen as a 0-1 programming problem, is NP-hard (Schuetz et al., 2022; Gasse et al., 2022). Therefore, finding efficient methods to relax the indicator matrix for optimization important.
Ng et al. (2001) relaxed the indicator matrix to the Steifel manifold, F E {X X = l}. where I the identity matrix. This approach further developed spectral graph theory and led the formulation of classic algorithms such as spectral clustering (Balestriero & LeCun, 2022; Macgregor & Sun, 2022). However, optimizing over the Steifel manifold always requires O e) operations (Wen & Yin, 2013), making it challenging to scale for large datasets, and it can only provide an optimal solution for problems of the form tr tr(F²2 LF), while in clustering, the resulting F still needs post-processing through methods like K-means (Li et al., 2015; Mondal al., 2021). An alternative relaxation is to make onto the single stochastic manifold, F E (X] X1, = In, X > 0} (Sun et al., 2015), which gave rise to well- algorithms like Fuzzy K-means (Ferraro, 2024; Borlea et al., 2021). However, this approach has the drawback of not considering the total number of samples per class, which can lead to empty clusters imbalanced class distributions (Ikotun et al., 2023; Hu et al., 2023). The most recent method is to relax the indicator matrix onto the double stochastic manifold, i.e., F {X X1. 1,, X²1, r, X 0} (Fettal et al., 2024; Yuan et al., 2024c). However, this approach also has significant drawbacks. The double stochastic manifold imposse overly striet requirements on the columns of F, as necessitates knowing the true distribution of
"Corresponding author.

---

Published as a conference paper at ICLR 2026
each class in the dataset as prior, which is nearly impossible for unknown datasets. Additionally, optimization over the double stochastic manifold is extremely challenging. still requiring O time (Douik & Hassibi, 2019; 2018), making almost infeasible for large-scale datasets.
To solve above questions, we propose a new relaxation method, where F E {X X1, = 1,, < 1",, < X 0). In this approach, the constraints the column sums are relaxed to lie within specified range. This allows us to flexibly incorporate as much prior knowledge as possible into the model. When there is more prior knowledge, we can choose a tighter (1, u) interval. Conversely, we can make more relaxed. Specifically, when the column sums and the true distribution are known, we can set u and to the true distribution (In fact, this does not lead to the absence of solutions. for further discussion, you see Appendix G). When no prior knowledge is available, we can set I < O and u >> n, which means our relaxation a generalization of both the single stochastic manifold and the double stochastic manifold, offering more adaptable framework.
We prove that the set of relaxed indicator matrices forms a manifold, which we call the Relaxed Indicator Matrix Manifold (RIM manifold). Based on Riemannian geometry (Boumal, 2023; Fei et al., 2025), we have developed a Riemannian optimization toolbox (Boumal et al., 2014; Townsend et al., 2016) for running optimization on the RIM manifold. In particular, we provide three distinet Retraction methods, including one that allows for fast computation of geodesics (Nguyen, 2022; Jordan et al., 2022), enabling algorithm to efficiently operate along the geodesic. Furthermore, we demonstrate that our algorithm, compared to existing Riemannian optimization methods on the double stochastic manifold, reduces the time complexity from O(n)) ) to (n). Furthermore, we have developed various Riemannian optimization algorithms that run on the RIM manifold.
We designed series of large-scale experiments with millions of optimization variables to validate our algorithm. These experiments include comparisons with state-of-the-art optimization algorithms on both convex and non-convex problems like image denoising (Takemoto et al., 2022; Zhou et al., 2024). In particular, we applied the Ratio Cut model (Veldt, 2023; Hagen & Kahng, 1992) to the RIM manifold. When u, our algorithm is 70-200 times faster than those based on the double stochastic manifold for large-scale problems with millions of variables, and it achieves lower loss results. In general, the algorithms on the RIM manifold outperform the latest optimization algorithms both loss function values and time. Additionally, the Ratio Cut clustering metric on the RIM manifold exceeds that the latest clustering algorithms.
Overall, our contributions include:
We propose a novel relaxation method for the indicator matrix, which allows for the full utilization of varying levels of prior information from the dataset, and we proved that the relaxed matrix forms manifold.
We develope a Riemannian optimization toolbox for manifolds, providing three Retraction algorithms, including a fast method for obtaining geodesics on the RIM manifold. We also demonstrated that the RIM manifold can replace methods on the double stochastic manifold, reducing the time complexity from Onn²) to (n).
We conducte lots of experiments with millions of variables, demonstrating the speed and efficiency of our algorithm. Our method outperforms the double stochastic manifold by 70-200 times in large-scale experiments, yielding better results and shorter time on various problems compared to latest optimization methods. We apply the RIM manifold to Ratio Cut and achieve superior clustering performance compared to the state-of tfate-of-the-art methods.
2 PRELIMINARIES
The Preliminaries section consists of four parts: an introduction to the notations, a brief overview of Riemannian optimization, and introduction to the single stochastic manifold, double stochastic manifold, and Steifel manifold, as well as machine learning methods on these manifolds. All the notations used in this paper follows the standard conventions of Riemannian optimization, and important symbols are introduced in the main text. Due to space limitations, the Preliminaries can be found in Appendix B.
2

---

Published as a conference paper at ICLR 2026
3 RIEMANNIAN TOOLBOX
3.1 DEFINITION OF THE RELAXED INDICATOR MATRIX MANIFOLD
The optimization of indicator matrix F Ind"xe where typically >> is NP-hard optimization problem. Three relaxation methods have already been introduced. The Steifel manifold E {X X X = } always requires O(n²c) time complexity (Wen & Yin, 2013) and can only yield an analytical optimal solution in the form of tr(F" LF),, while in clustering, the resulting F still needs post-pree post-procsssing through methods like K-means. The single stochastic manifold
{X X1. In. X 0} does not impose any constraints on the column sums of F, which may lead to empty or imbalanced classes and cannot incorporate column sum information into the model. The double stochastic manifold F E {X] X1, n, X 1, time complexity of (n²), and the constraints on the column sums are too strict, often making impossible to obtain the sum of the column. Therefore, we propose new relaxation method:
X > 0), on the other hand, has
F {X XLc = Ln, < 1 u, > 0}
(1)
Introducing / and u allows us incorporate much information as possible into the model. Ad-
ditionally, when < 0 and > n, our relaxation reduces to {X Xle In, X > Oj. When r, our relaxation becomes {X X1. In, 1 r, 0}. Thus, our relaxation generalizes the previously mentioned approaches. Importantly, our relaxation forms an embedded submanifold of the Euclidean space.
Lemma 1. Our relaxed indicator matrix set == {X XI, = -,, X1,, < u, X > 0f forms
an embedded submanifold of the Euclidean space, with dim M (n 1)c. We refer to it as the Relaxed Indicator Mairix Manifold. The proof is included in A.
3.2 RIEMANNIAN OPTIMIZATION TOOLBOX FOR THE RIM MANIFOLD
In this section, we will establish an optimized Riemannian toolbox for the RIM manifold. To transform the embedded submanifold (Zhang et al., 2024; Lee & Lee, 2012) M into a Riemannian submanifold (Lee, 2018; Gulbahar, 2021), it is necessary to equip MM with an inner product (.-)xx- Mishra et al. (2021) adopt the Fisher information (Ly et al., 2017; Rissanen, 1996) metric for manifolds. However, an alternative approach is to directly restrict the Euclidean inner product onto the manifold. The reason for doing is seen in F. This restriction allows for straightforward derivation of the Riemannian gradient (Huang & Wei, 2022) from the Euclidean gradient and the method lies in enabling an intuitive and convenient Retraction mapping.
Lemma 2. By restricting the Euclidean inner product (U, = En Σ Vig onto the RIM
manifold M, the tangent space of M at X is given by Tx (U Ul 0}. For any function H. if Euclidean gradient Grad H(F).. the Riemannian gradient grad H(F) expressed as following. The proof is included in A.2
grad H(FF- = Grad H(FF) 1 Grad H(F)).1T
(2)
To further obtain second-order information of a function, it is necessary to equip the manifold M with Riemannian connection (Epstein, 1975). We select the unique connection that ensures the Riemannian Hessian hess H is symmetric and compatible with the inner product as the Riemannian connection. The following theorem formalizes this:
Lemma 3. For the manifold M, there exists a unique connection that is compatible with the inner product and ensures that the Riemannian Hessian mapping is self-adjoint. This connection is given by following. V is the Riemannian connection in Euclidean space. The proof is included in A.3
The Riemannian Hessian mapping can be directly derived from the above Riemannian connection. Lemma 4. For the manifold M equipped with the connection VvU, the Riemannian Hessian mapping satisfies following. Hess H is the Riemannian Hessian Euclidean space. The proof is included in A.4
1 hess H[V] = Hess H[V1 Hess H[V]1.17
(4)
3

---

Published as a conference paper at ICLR 2026
A Retraction (Boumal, 2023; Absil al., 2008) is mapping Rx (tV) that maps from the tangent space of M at X to the manifold M, i.c., Rx - 4 A Retraction is used to generate curve y(t) Rx (tV), starting X and moving in the initial direction given by V. allowing X to move along the manifold. Specifically, Rx (tV) should satisfy Rx (0) X and d Rx((t)) --: V. ăt
D dt (0)) = O, then 7(1) forms geodesic, where dt represents the Levi-Civita derivative (Berz, 1996). Geodesics provide better convergence guarantees for optimization algorithms manifolds (Vishnoi, 2018). The following theorem presents a method for obtaining geodesics.
t=u
Theorem 1. Let Rx (tV) argmin
(X + tV X E M. Then
(5) where v(t), (t)) (t) are Lagrange multipliers. Moreover, there exists 8 such that for t E (0. 6), - (,w(t)T + Imp(t)" O, and the Retraction satisfies the following. Where D at
argmin FEM (X + + max(0, X + tV v())1? - 1, (t) 1np (())((t)
denotes the Levi-Civita derivative.
Rx(0) Xx d
Rx =V,
D at R =0
di
1=
Thus, Rx (IV) is a geodesic. The proof is included in A.5.
The essence of solving the Retraction is to compute the orthogonal projection geminngeMll F (X II which be addressed from two perspectives: the primal problem and the dual problem.
Theorem 2. MM 21 2gg n llg where Q2 (X X X1, 1n). S2 {X X l}, and Q2 = {X 1 < u}. The primal problem can be solved using the Dykstras (Tibshirani, 2017; Boyle & Dykstra, 1986) algorithm by iteratively projecting onto D2 f2 and Zz. Specifically:
Projo, (X) = (Xij + 4 where is determined by Projo, (X)le = 1m-
Projaz (X) and Proja, (X) are defined similarly. For example,
[X), L1 IX),,, +X, if ( S 44
if ( ر4
Projaz (x))
where X is the j-th column of X, and 1j is the j-th element of the column vector l. The proof is included in A.6.
Please note that Dykstras algorithm is projection based based retraction algorithm (Absil & Malick, 2012). which in fact seeks the best point in the set {X] X1, = 1,, S X",, S u, X 2 0}. This design allows the algorithm to remain effective even when 1 = A detailed discussion of this can be found in Appendix G.
Another approach is the dual gradient ascent method. We have proven the following theorem. Theorem 3. Solving the primal problem is equivalent to solving the following dual problem:
1 _ max X + vv 1,w"" 1np")) (v, In) (w, u) + (p.1)
max C # 200.20
(8)
where v, w. and p are Lagrange multipliers. The partial derivatives of L with respect to v, w, and p are known, and gradient ascent can be used solving v, w, and p. Finally, Rx(tV) (tV) can be obtained using max(0, X tV v1 Inw"" + +Inp²') The partial derivatives are following. The proof included in A.7.
dv m0 X +V v1 1,, + 1np" 11 l
DL
dw m0 tV v1 1, + 1mp 7) F
ac max(0, X + tV - -' - 1n," T + Lap 7 +1
0p
Additionally, propose Retraction method based variant of the Sinkhorn algorithm (Xie et al., 2025; Cuturi, 2013). This approach also attempts to map a matrix onto the RIM manifold using two diagonal matrices. The following theorem illustrates this property. However, it is equivalent to solving an optimal transport problem with an entropy regularization parameter, whose choice may not be well justified.
4

---

Published as a conference paper at ICLR 2026
Table 1: Time complexity comparison(n (n > c).
Operation
Additions
RIM Manifold
Multiplication
Total
Additions
Deuhly Stochastic Manifold
Multiplications
Total
Speedup factor
Riemannian Gradient
D(nc]
O(n)
Oia)
Oin")
Oja")
Ola²)
Retraction
Oquc)
Oinc)
Oluc)
O(nc)
Oine)
Oil)
Riemannian
Olyc)
((n)
O(n)
Olr²)
Oin")
O(a²)
O(w²)
Theorem 4. The Sinkhornbased Retraction is defined as (XV(VV) S(X @ exp(IV @ x)) diag(p") (XX o exp(lV @ X)) diag(a' @ w=") where p',q", vector into diagonal matrix. The vectors p" o w are obtained by iteratively updating the following equations:
(10)
, w
are vectors, exp(-) denotes element-wise exponentiation, and diag(-) converts
pp(+1) = =I, 2 ((X 3 expltV @ x)) (g(h)) o w(())),
(b-1)) = mmxx 110 (X @ exp(tV @ x))² plot1) @ will)) 1e): g
(11)
w((+1) min u (X exp(tV ex) pltt e qlatt)) 1o).
This iterative procedure ensures the mapping onto the RIM manifold. The solution R& (IV) diag(p") )X @exp(tV (x) diag(a" @w") is equivalent to solving the dual-bound optimal transport problem (12) with an entropy regularization parameter of I. The proof is included in A.8.
Rx (LV) argming ( log(X @ @exp(IV @ xI)) + 8 6-6=
n
(Feg (og()) Fig) (12)
Based on the Riemannian toolbox for the RIM manifold, we have developed Riemannian Gradient Descent (RIMRGD), Riemannian Conjugate Gradient (RIMRCG), and Riemannian Trust-Region (RIMRTR) methods on the RIM manifold. The algorithmic procedures are provided Appendix C.
3.3 COMPARISON ANALYSIS OF TIME COMPLEXITY
When 1, the RIM manifold reduces to the doubly stochastic manifold and provides a fast way for solving problems on the doubly stochastic constraint. Existing optimization methods on the doubly stochastic manifold are extremely time-consuming. This section provides comparative analysis of the time complexity between the RIM manifold and the doubly stochastic manifold.
First, discuss the Riemannian gradient. The computation of the Riemannian gradient the RIM manifold is given by grad ((F) = Grad H(F) - 1/2 Grad H(F)1,17 Here, the term Grad H,F) involves summing each column, dividing by c. and then replicating it across c columns. This requires 2nc additions and divisions.
For the doubly stochastic manifold, the Riemannian gradient (n = c):
grad H(F) = - (a17 + 1,117 - 1,aa F) ① F, a= ((-FFE)) (7
(13)
- (
Grad P((F') O F.
The term FFT Raxx, and computing its pseudo-inverss (I- FFTY requires least n additions or multiplications. Further computing the Riemannian gradient involves at least n² operations. When n * e, we need to solve linear system of (n + e) dimensions still takes ( time (where n >>>>>> c).
For the Retraction operation, the time complexity is O(nc), which scales linearly with the number of variables. For the computation of the Riemannian Hessian, the RIM manifold also requires only (nc) additions and Occ) multiplications. In contrast, the Hessian mapping the doubly stochastic manifold has a highly complex expression (186), requiring at least (n additions and multiplications.
We summarize the time complexity in Table 1, including the complexity of each operation and the speedup factor. We will conduct extensive experiments to verify the acceleration effect.
4 RIM MANIFOLD FOR GRAPH CUT
In this section, we apply the RIM manifold to graph cut problems, using Max Cut (Shinde et al., 2021; Wang et al., 2022) and Ratio Cut (Chen al., 2022h; Nie et al., 2024) as examples. Max Cut
5

---

Published as a conference paper at ICLR 2026
and Ratio Cut are both well-known graph partitioning algorithms, and their loss functions given by Hm((P) -tr(F" for the Max Cut, and (F(F)= = tr tr(F"LE(F"" -)- for the Ratio Cut. S is the similarity matrix, and L is the Laplacian matrix (Nie et al., 2016; 2014). The constraint is F E Ind"xe, and we relax this constraint on the RIM manifold.
First, the Euclidean gradient of -tr(E" SF) is Grad(-tr(F" Riemannian gradient grad Hm (F) = SF SF1,1 According to Theorem 4, the Riemannian Hessian expression is hess H. [V] = Hess Hm [V] [m that:
= SF,, and its corresponding
- Hess Hm [V]122 Moreover, we know
Grad Hm (F + (V) Grad (m
- S(F + (V) SF = -V
Hess 7m V
lim z-0
lim
(14)
t
Therefore, show that hess
(r((F" can be represented as following:
= -V --VILL
hess
(15)
Now we apply the RIM manifold to the Ratio Cut problem. Ratio Cut is an important graph partitioning method with the objective function (((FTFF((F) FF)(FF))))) subject to F E Indexe The relaxed optimization problem is formulated as:
min tr(F" LF(F" F)""), M
(X X1. L., < x u, X 0
(16)
FEM
The following theorem provides the expressions for the Euclidean gradient and the Euclidean Hessian map of the Ratio Cut.
(EF(FFF)-))) Theorem 5. The loss function for the Ratio Cut is given by H,(F) = tr(FF LF(FT F)-2), Then, the Euclidean gradient of the loss function with respect to F is following. The proof is included in A.9.
GradH (F) = - (( F)))
p) F) -
(17)
Given the substitutions (F² F)-1 J and F LF K, the Euclidean Hessian map for the loss function is:
HessHo = 2 - F F"V - VJKJ + FJ(V" F F F'VJJKJ
(18)
- -FJ(V" F +F LV)J (OJJ+FJJJY" F"V)J)
(19)
The above theorem provides the Euclidean gradient of Ratio Cut. Although computing (F" (FTF)-1 requires inversion, where F E Re the inversion complexity is only O(c²) and c < n. Next, we will perform graph optimization the RIM manifold, comparing the loss results and runtime with various state-of-the-art algorithms, as well as evaluating the effectiveness of graph cut for clustering.
We further prove that the graph problem defined on the RIM manifold is always Lipschitz continuous. More specifically, we have the following theorem.
Theorem 6. For any graph cut problem expressed as H(FF) tr((FT LF) is any symmetric matrix, the Fuclidean gradient GradH (F) is bounded, and satisfies:
where W
[GradH(F))l 22 [ILIoVV |W
mmin (W)- -2
(20)
a²
where anin (W) is the smallest singular value of the matrix W and 11-11 denotes the spectral norm. This implies that H(F) is Lipschitz continuous.
In addition, we provide convergence theorems for graph cut optimization on the RIM manifold using Riemannian optimization. Specifically, in Proof A. we show that the graph cut problem is Lipschitz smooth and provide an analysis of its convergence rate. We also analyze the relationship between finding the optimal solution via the doubly stochastic manifold and the RIM manifold under the assumptions of strong convexity and L-smoothness. The proof can be found in Appendix A.
5 EXPERIMENTS
In this section, we will conduct extensive experiments to evaluate the performance of Riemannian optimization on the RIM manifold and address several key questions of interest.

---

Published as a conference paper at ICLR 2026
Question 1: For the RIM manifold, this paper proposes three different Retraction methods. Which method is the most efficient? Which Retraction is recommended for use?
When u. does the Riemannian optimization algorithm the RIM manifold outperform the Riemannian optimization algorithm on the doubly stochastic manifold in terms of effectiveness and speed?
Question 2:
Question 3: For non-convex optimization problems, we evaluate whether optimization on the RIM manifold is faster or more effective compared to other state-of toffthe-artt methods? As examples, we consider a classic non-convex graph cut problem Ratio Cut.
Question 4: When relaxing the graph cut problem onto the RIM manifold (followed by discretization), can common clustering metrics(ACC,NMII.AR achieve better values?
5.1 EXPERIMENTAL SETUPS
5.1.1 EXPERIMENT SETUP
To determine which of the three Retraction methods is more efficient, we randomly select a large number of matrices V M,, i.e., generate large number of tangent vectors, and set t = 1. Then, we apply the three Retraction methods to generate points on the RIM manifold M. To ensure the experiment's 's validity, from (500, 1000, 3000, 5000, 7000, 10000} and e takes values from (5. 10, 50, 100, 500, 1000). The lower and upper bounds are set as = 0.9 [ and v 1.1 [ J. respectively, as well as ( /= = = We then calculate the computation time for the three Retraction methods and compare them. For large-scale problems, we recommend using the faster Retraction method. the computation times nearly identical, recommend using the norm- norm-based Retraction, as it yields geodesics with better properties.
vary the matrix dimensions V E Ra where takes values
5.1.2 EXPERIMENT 2 SETUP
To answer the second question, we need to compare Riemannian optimization methods on the RIM manifold with optimization methods the doubly stochastic manifold under the condition / = To this end, we design two optimization problems, including both convex and non-convex cases.
The first problem is norm approximation problem. Specifically, E "x with sizes n E (5000, 7000, 10000), E {5, 10, 20, 50, 70, 100) and solve the following optimization problem. We compare the runtime and loss function values of the two manifolds.
randomly generate matrix
min FEM
|F
M (X X1. Ale, A 0}
(21)
The second problem is image denoising task based on the classical total variation (TV) regulariza- tion model. The RIM-TV model is given by
minpeM 1 All
Eii (F - Fi, + [F+1., Figl) M- (X X 0, X1. = =Al,, X = 1"
(22)
Here, & is the total variation (TV) regularization coefficient, A is the original image obtained from the dataset, and Á is the noisy image generated by adding Gaussian white noise to A. The image A; i0, 1). E is chosen from the set [0.3,0.7), and the variance of the added Gaussian noise is chosen from the set {0.3, 0.5, 0.9]. We will compare the speed and objective function values of the algorithm when running on RIM manifold versus the doubly stochastic manifold. More experimental details can be found in Appendix D. When comparing the RIM manifold with the doubly stochastic manifold, the only difference lies in the manifold toolbox used, all other components, including the line search method and stopping criteria, kept exactly the same.
5.1.3 EXPERIMENT 3 SETUP
To answer the third question, we apply the RIM manifold to Ratio Cut and conduct experiments on real datasets (as shown in Appendix D.3.4). The values and u are set as / = u = and / = 0.9 _ u = 1.1 n respectively. For = u = g. we compare seven algorithms: Riemannian Gradient Descent (RIMRGD), Riemannian Conjugate Gradient (RIMRCG), Riemannian

---

Published as a conference paper at ICLR 2026
Table 2: Table of Execution Time when / = u for Different Retraction Algorithms(
Row&Col
500
1000
3030
Daal
2000
10000
1000
2000
Sinklon
1200
3000
Dykstzas
5200
0,012
0.024
0.106
0.140
0.00)
0.043
0.007
0.011
10
0.022
0.036
0.075
0.112
0.641
0.188
9,802
0.005
0.036
0.087
0.166
0.343
0.036
0.985
0.814
0.833
0.831
0:843
50
0.074
0005
0.791
1:307
1,866
3.766
9,802
0,805
01641
0.175
0.278
0.674
0036
0
0.823
0.89
0.851
0.874
100
0.012
0.174
1.997
2,922
3
5.710
9.003
8,805
0.006
0.140
0.258
0.580
0.966
0,010
0.831
0,800
0.872
0.306
509
0.994
6,813
0.830
0.593
01629
0006
0.1%
0.890
0.168
1000
0.102
0.178
17.36
25.56
40.55
36.56
0/034
0,082
0.446
3.14
0,867
0.319
0.364
Trust Region (RIMRTR), Frank-Wolfe Algorithm (FWA) (Jaggi, 2013; Weber & Sra, 2023; Yurtsever
& Sra, 2022), Projected Gradient Descent (PGD) (Shen & Chen, 2023; Chen & Wainwright, 2015), Riemannian Gradient Descent on the Double Stochastic Manifold (DSRGD) (Tripuraneni et al., 2018), and Riemannian Conjugate Gradient on the Double Stochastic Manifold (DSRCG) (Sato, 2022). For 1 0.9 [ and PGD. The optimization results of these algorithms are then compared. More experimental details can be found in Appendix D.2.
2 and 1.1 ]. only compare RIMRGD, RIMRCG, RIMRTR, FWA,
5.1.4 EXPERIMENT 4 SETUP
To answer the fourth question, we compare the Ratio Cut algorithm on the RIM manifold with ten clustering algorithms. We again choose 8 real datasets with different types, including images, tables, waveforms, etc. (as shown in Appendix D.3.4), and conduct large-scale validation using 10 comparison algorithms (listed in D.3). We evaluate the clustering performance using three metrics: clustering accuracy (ACC) (Yuan et al.. 2024a;b), normalized mutual information (NMI) (Ren et al., 2024), and adjusted Rand index (ART) (Ronen al., 2022). For the similarity matrix, k-nearest neighbor (k-NN) (Li et al., 2024b; Zhu et al., 2022) Gaussian kernel function (Wang al., 2009; Chen et al., 2021) and construct the Gaussian kernel function using the mean Euclidean distance. For the parameter k. each comparison algorithm is tested by searching for the best value of within the range [8, 10, 12, 14, 16). More experimental details can be found in Appendix D.3.
the
5.2 EXPERIMENTAL RESULTS
5.2.1 RESULT OF EXPERIMENTAL.
The data for Experiment I when / = is presented in Table 2. The horizontal axis indicates the methods used, while the vertical axis represents the number of columns, and the horizontal axis represents the number of rows of the experimental matrix. The table entries represent the time required for Retraction, measured in seconds. The fastest method is highlighted in red. As observed, when the matrix is small, the Sinkhorn method is faster. However, the matrix size increases, the Dykstras method shows significant advantages and produces the geodesic. Therefore, we recommend using the Dykstras method obtain the Retraction curve. More data can be found in Appendix E.
5.2.2 RESULT OF EXPERIMENTAL 2
Table shows the time and final loss required by the Riemannian Trust Region method to solve convex optimization problems of different scales. It can be seen that, for problems of varying sizes, the RIMTRT significantly outperforms the DSTRT in both time consumption and final loss. Therefore, we have highlighted the RIM manifold results in red. Data for the Riemannian Gradient Descent and Riemannian Conjugate Gradient methods can be found in Table 13 and Table 14.
For the second part of the experiment, Figure 1 shows the comparison of denoising results using the TV algorithm on the RIM manifold and the doubly stochastic manifold with a noise level of 0.3. In this case, { = 0.3. On the RIM manifold, the running time was 29.77s, and the loss value decreased to 1.05e5, while on the doubly stochastic manifold, the time was 85.33s, and the loss value was 1.17e5. By observing the images, it is evident that the image obtained using the doubly stochastic manifold has noticeable noise when zoomed in, while the image on the RIM manifold is smoother. Additional data and images can be found in Figure
In addition, we evaluated the denoising results obtained from the RIM manifold and the doubly stochastic manifold using multiple metrics ( (Valsesia et al., 2020). The resulting table is shown Table 4. It can be observed that the denoising performance of the RIM manifold consistently surpasses that of the doubly stochastic manifold. In addition, it is important note that the relative values rather
8

---

Published as a conference paper at ICLR 2026
Table Cost and Time on the RIM Manifold and Doubly Stochastic (RT)))l
Row&Col
RIM Manifold
Doubly Stochastic Manifald
Cost
Time
Cost
Time
10000
3.09E-23
8.288-20
2,091:220
0.265
0.355
0.516
4,38E-11
3,96E-10
2.89E-10
9.530
12.85
31.97
10
1. 1.908-200
9.588-20
3 3.808200
0.2283
0.464
0.660
1.918-10
3.668-10
4.125-10
16.25
17.04
32.72
1.02E-19
1,228-23
8.295-19
0.366
0.562
9,49B-10
6.04B-10
1.17E-499
7.66E-20
2.136.18
2.206-20
3.08E-09
1.99E.09
70
1 1.858-20
2.84x-18
7.498-19
0.791
0.983
1.352
2.59E-499
1.55E-499
3.18F-09
70.47
121.0
77.28
100
1.31E-19
5.048-20
1,268-17
0,950
1.324
1.721
1,788-49
2.18E.499
2.83E-099
91.40
121.3
241.4
(a) Origin
(b) Noisy Image
(c) RIM Result
(d) DS Result
Figure 1: Image Denoising Results, Noise Coefficient 0.3, 5 = 0.3.
than the absolute values- metries such LPIPS should be the focus, since TV regularization is not a state-of-the-an denoising method. Our goal here is to compare RIM against DSM.
5.2.3 RESULT EXPERIMENTAL
When u, the time and loss for the seven comparison algorithms are presented Table We have marked the algorithm names on the RIM manifold in blue, the shortest time in red, and the lowest loss bright red. can be observed that the optimization algorithms on RIM manifold achieved most of the top positions. Figure 2 shows the loss decrease curves for some datasets. More results can be found in Appendix E.3.
5.2.4 RESULT EXPERIMENTAL
For Experiment Table records the performance of comparison algorithms across real- world datasets based on clustering accuracy (ACC), normalized mutual information (NMI), and adjusted Rand index (ARI). Our algorithm marked in blue, and the best-performing algorithm marked in red. It can be observed that performing Ratio Cut on the RIM manifold leads to superior results compared to the most advanced algorithms. More results can be found in Appendix E.4.
LIMITATIONS
We acknowledge that our study still has several limitations that warrant further investigation. First, the RIM manifold a relaxation of the indicator matrix, and for an NP-hard non-convex indicator-matrix optimization problem, relaxations such as the single stochastic, doubly stochastic, and Stiefel manifolds various aspects. It remains challenging provide precise bound between the optimal solution of the relaxed problem and that of the original NP-hard problem under general conditions. Second, the RIM manifold its retraction. When strict adherence to open set is required,
analysis only demonstrate that the RIM manifold outperforms existing
employs projection onto closed a small correction term e. Finally, additional experimental results could be included in future work to further demonstrate the superiority the RIM manifold.
(a) USPS20
(b) JAFFE
(c) MnistData05
(d) waveform21
Figure 2: Comparison of Loss Decrease for Optimization Algorithms on Real Datasets.

---

Published as a conference paper at ICLR 2026
Table Comparison of Denoising Performance between the RIM Manifold and the DSM
(noise, E Metric
0..3.0.3
(0.3. 0.73
(0.5.
(0.5.
0.7)
(0,90.33
(0.9,
DSM
RIM
DSM
DSM
MSE
0.012
0.015
0.020
0.027
0.022
0.026
0.023
0.032
0.085
0.106
0.044
0.051
PSNR
19.27
9.751
SSIM
0.502
0.412
0.434
0.256
0.282
0.247
0.327
0.209
0.107
0.096
0.183
0.148
LPIP
0.561
0.671
0.719
0.929
0.742
0.824
0.775
0.968
0.968
1,020
0.803
1,020
Table 5: Time and Loss Different Optimization Algorithms on Ratio Cut when |
DSRGD
USRCG
FWA
PGD
KIMRGU
RIMRUG
RIMRTR
Time
Cost
Time
Cost
Time
Coo
Time
Clst
Time
Cos
Time
Cost
Timee
Clss
COIL20
3.97%
25.17
11.90
23.41
10.49
43.12
0.967
31.53
1.145
24.83
0.685
27.40
14.20
22.45
Digit
3.000
2.751
2.312
9.196
6.077
0953
7.058
0.983
1.119
13:75
1.089
2.725
29.56
2.775
2.221
1,070
516.3
Waszform22
2.300
MaisiDatal5
2.520
2.126
2.154
Table 6: Mean clustering performance of compared methods real-world datasets.
Metric
Method
COIL.0
Digit
JAFFE
MSRA25
PalmData25
USPS20
Waveform21
MnistData05
53.44
58,33
72.16
49.33
70.32
55.51
50.38
53.86
CDKM
52.47
65.82
80.85
59.63
76.05
57.68
50.36
54.24
Reut
78.14
74.62
84.51
56.84
87.03
57.83
51.93
62.80
Neut
78.88
76.71
83.76
56.23
86.76
59.20
51.93
61.14
Nystrom
51.56
72.08
75.77
52.85
76.81
62.55
51.49
55.91
BKNC
57.11
60.92
93.76
65.47
86.74
62.76
51.51
52.00
ACC
FCFC
43,94
$4.27
69.38
58.23
$6.98
$4.41
FSC
82.76
79.77
81.69
56,25
82.27
67.63
50.42
57.76
LSCR
65.67
78. 14
91.97
53.82
58.25
63.07
56.19
57.15
LSCK
62.28
78.04
84,98
$4.41
58.31
61.86
54.95
58.57
RIMRcut
79.72
82.53
96.71
56.64
90.85
70.28
74.80
65.55
KM
71.43
58,20
80,93
60.10
89,40
54,57
36.77
49.57
CDKM
71.16
63.64
87.48
63.83
91.94
55.92
36.77
49.23
Reut
86.18
75.28
90.11
71.64
95.41
63.84
37.06
63.11
Neut
86.32
76.78
89.87
71.50
95.26
64.46
37.06
63.22
Nystrom
66.11
70.13
82,53
57.77
93.09
59.00
36.95
48.53
BKNC
69.80
59.37
92.40
69.30
95.83
57.10
36.94
44.56
NMI
74.05
38.33
80.30
63.34
89.47
55.71
22.89
48.75
91.45
80,98
90.43
70.60
94.62
74.75
36.76
LSCR
74.67
75.07
93.13
68.06
81.84
62.36
33.37
52.82
LSCK
74.02
76.53
87.89
67.97
81.70
65.23
36.92
59.14
RIMRcut
85.63
80.05
96.24
71.76
96.50
69.08
42.14
59.35
KM
50.81
45.80
66.83
34.66
65.06
43.57
25.56
37.18
CDKM
48.11
52,74
76.36
37.70
71.73
45.59
25.56
36.79
73.73
81.70
46.35
84.76
$1.99
25.31
Neut
74.30
68.21
81.30
45.90
84.25
$2.72
25.31
50.51
Nystrom
45.96
59,50
69.85
38.07
76.23
50.01
25.03
38.21
BKNC
49.96
48.98
87.96
54.78
85.56
48.43
25.02
32.89
ART
FCFC
54.41
25.50
65.73
40.42
66.03
46.32
22,89
36.86
FSC
79.46
73.03
80.26
43.99
79.67
61.71
25.10
44.78
LSCR
67.21
43.31
48.70
$2.64
25.12
LSCK
54.59
68.70
77.37
42.18
48.58
52.54
26.47
46.48
RIMReut
73.98
75.01
93.32
46.82
88.49
56.06
42.89
52.87
CONCLUSION
This paper presents new relaxation for indicator matrices and proves that it forms a Riemannian manifold. We have constructed Riemannian toolbox for optimization on the RIM manifold. In particular, we introduce multiple methods for Retraction, including the geodesic. The paper demonstrates that optimization on the RIM manifold is useful for machine learining and it is fast method (nn) that replace the existing double stochastic manifold optimization with a time complexity of O(n²). Through large-scale experiments from multiple perspectives, have proven the effectiveness and speed of optimization the RIM manifold.
that operates quickly along
ACKNOWLEDGEMENT
This work supported by the National Natural Science Foundation of China under Grant 62576277.
10

---

Published as a conference paper at ICLR 2026
9 STATEMENT
For the reproducibility of this paper, we have submitted the complete anonymized code with fixed random seeds, detailed in Appendix In addition, large language models (LLMs) were only used for language polishing.
REFERENCES
Hervé Abdi and Lynne J Williams. Principal component analysis. Wiley interdisciplinary reviews: computational statistics, 2(4):433 459, 2010.
P-A Absil and Jérôme Malick. Projection-like retractions on matrix manifolds. SIAM Journal on Optimization, 22(1):135-153 158, 2012.
P-A Absil, Robert Mahony, and Rodolphe Sepulchre. Optimization algorithms on matrix manifolds. Princeton University Press, 2008.
Randall Balestriero and Yann LeCun. Contrastive and non-contrastive self-supervised learning recover global and local spectral embedding methods. Advances in Neural Information Processing Systems, 15:26671-26688 2022.
Yichen Bao, Han Lu, and Quanxue Gao. Fuzzy k-means clustering without cluster centroids. arXiv preprint M.XXv:2444400944: 2024.
Jean-David Benamou, Guillaume Carlier, Marco Cuturi, Luca Nenna, and Gabriel Peyré, Iterative bregman projections for regularized transportation problems, 2014. URL https://arxa
h/g//bb//4422 5154.
Martin Berz. Calculus and numerics on levi-civita fields. Computational Differentiation: Techniques, Applications, and Tools, (89):19-37, 1996.
James Bezdek, Robert Gunderson, Robert Ehrlich, and Tom Meloy. On the extension of fuzzy k-means algorithms for detection linear clusters. In 1978 IEEE Conference on Decision and Control including the 17th Symposium on Adaptive Processes, pp. 1438 1443. IEEE, 1979.
Ioan-Daniel Borlea, Radu-Emil Precup, Alexandra-Binnn Borlea, and Daniel lercan. A unified form of fuzzy c-means and k-means algorithms and its partitional implementation. Knowledge-Bass Systems, 104:10 2021.
Nicolas Boumal. Optimization and estimation on manifolds. 2014.
Nicolas Boumal. An introduction to optimization on smooth manifolds. Cambridge University Press, 2023.
Nicolas Boumal, Bamdev Mishra, P-A Absil, and Rodolphe Sepulchre. Manopt, matlab toolbox for optimization on manifolds. The Journal of Machine Learning Research, 15(1):1455-1455 2014.
Stephen Boyd and Lieven Vandenberghe. Convex optimization. Cambridge university press, 2004.
James P Boyle and Richard I. Dykstra. A method for finding projections onto the intersection of convex sets hilbert spaces. Advances in Order Restricted Statistical Inference: Proceedings of the Symposium Order Restricted Statistical Inference held in Iowa City, Iowa, September 11-13, 1985, pp. 28-47. Springer, 1986.
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 (SampTA), pp. 73-77. IEEE, 2017.
Xiaojun Chang, Feiping Nie, Zhigang Ma, and Yi Yang. Balanced k-means and min-cut clustering, 2014. URL pattps://arxx
woo.o///baa.............. abs /111.. 6235.
Lactitia Chapel, Mokhtar Z. Alaya, and Gilles Gasso. Partial optimal transport with applications on
positive-unlabeled learning, 2020. URL /ttetpss//axxi..rrrrr/ httrggabs//0000 08276 6.
11

---

Published as a conference paper at ICLR 2026
Huimin Chen, Qianrong Zhang, Rong Wang, Feiping Nie, and Xuelong Li. A general clustering framework based on a novel balance regularizer. Signal Processing, 198:10,57, 2022a.
Yen-Yen Chen, Yangqiu Song, Hongjie Bai, Chih-Jen Lin, and Edward Y. Chang, Parallel spectral clustering distributed systems. IEEE Trans. Puttern Anal. Mach. Intel., 33(3):668-588 2011.
Xiaojun Chen, Zhicong Xiao, Feiping Nie, and Joshua Zhexue Huang, Finc: An efficient and effective optimization method for normalized cut. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2022b.
Xinlei Chen and Deng Cai. Large scale spectral clustering with landmark-based representation. In Proceedings of the Twenty- AAAI Conference on Artificial Intelligence, AAAI 2011, San Francisco, California, USA, August 7-11, 2011. AAAI Press, 2011.
Yifan Chen, Qi Zeng, Heng Ji, and Yun Yang. Skyformer: Remodel self-attention with gaussian kernel and nystrom method. Advances in Neural Information Processing Systems, 34:2122-2133 2021.
Yudong Chen and Martin J Wainwright. Fast low-rank estimation by projected gradient descent:
General statistical and algorithmic guarantees. arXiv preprint arXiv: 03025, 2015.
Marco Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. Advances in neural information processing systems, 26, 2013.
Zhiyuan Dang, Cheng Deng, Xu Yang. Kun Wei, and Heng Huang. Nearest neighbor matching for deep clustering. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pp. 13693-13700 2021.
Ahmed Douik and Babak Hassibi. A riemannian approach for graph-based clustering by doubly stochastic matrices. In 2018 IEEE Statistical Signal Processing Workshop (SSP), pp. 806-810. IEEE, 2018.
Ahmed Douik and Babak Hassibi. Manifold optimization over the set of doubly stochastic matrices: A second-order geometry. IEEE Transactions on Signal Processing, 67(22):5761-5774, 2019.
Alan Edelman, Tomás A Arias, and Steven T Smith. The geometry of algorithms with orthogonality constraints. SIAM journal on Matrix Analysis and Applications, 20(2):003-33,, 1998.
DBAd Epstein. Natural tensors riemannian manifolds. Journal of Differential Geometry, 10(4): 631-645, 1975.
Jicong Fan, Yiheng Tu, Zhao Zhang, Mingbo Zhao, and Haijun Zhang. A simple approach to automated spectral clustering. Advances in Neural Information Processing Systems, 35:9907-9921. 2022.
Yanhong Fei, Yingjie Liu, Chentao Jia, Zhengyu Li, Xian Wei, and Mingsong Chen. A survey of geometric optimization for deep learning: from euclidean space to riemannian manifold. ACM Computing Surveys, 57(5):1-37, 2025.
Maria Brigida Ferraro. Fuzzy k-means: history and applications. Econometrics and Statistics, 30: 110-123, 2024.
Chakib Fettal, lazhar labiod, and Mohamed Nadif. Graph cuts with arbitrary size constraints through optimal transport. Transactions on Machine Learning Research, 2024. ISSN 2835-8856. URL https:
htttps://openreve wtwp/eeeeevvieeeeeeeeeeee/eeeeee ew net / forum? d=UG7 rt rsuaT.
Kyle Fox, Debmalya Panigrahi, and Fred Zhang. Minimum cut and minimum k-cut in hypergraphs via branching contractions. ACM Transactions on Algorithms, 19(2):1-25 2023.
Maxime Gasse, Simon Bowly, Quentin Cappart, Jonas Charfreitag, Laurent Charlin, Didier Chételat, Antonia Chmiela, Justin Dumouchelle, Ambros Gleixner, Aleksandr M Kazachkov, al. The machine learning for combinatorial optimization competition (ml4co): Results and insights. In NeurlPS 202/ competitions and demonstrations track, pp. 220-231. PMLR, 2022.
12

---

Published as a conference paper at ICLR 2026
Mehmet Gulbahar. Qualar curvatures of pseudo riemannian manifolds and pseudo riemannian submanifolds. AIMS Mathematics, 6(2):1366-1377 2021.
Lars Hagen and Andrew Kahng. New spectral methods for ratio cut partitioning and clustering. IEEE transactions computer-aided design of integrated circuits and systems, 1109):107--0085, 1992.
Wei He, Shangzhi Zhang, Chun-Guang Li, Xianbiao Qi, Rong Xiao, and Jun Guo. Neural normalized cut: A differential and generalizable approach for spectral clustering, Pattern Recognition, 164: 111545, 2025.
Haize Hu, Jianxun Liu, Xiangping Zhang, and Mengge Fang, An effective and adaptable k-means algorithm for big data cluster analysis. Pattern Recognition, 139:109404, 2023.
Wen Huang and Ke Wei. Riemannian proximal gradient methods. Mathematical Programming. 194 (1):371-413, 2022.
Mia Hubert, Peter Rousseeuw, and Karlien Vanden Branden. Robpca: new approach to robust principal component analysis. Technometrics, 47(11)644-7,, 2005.
Abiodun M Ikotun, Absalom E Ezugwu, Laith Abualigah, Belal Abuhaija, and Jia Heming. K-means clustering algorithms: A comprehensive review, variants analysis, and advances in the era of big data. Information Sciences, 622:178-2161 2023.
Martin Jaggi. Revisiting foank-wolfe: Projection-free sparse convex optimization. In International
conference on machine learning, pp. 427-435, PMLR, 2013.
Bo Jiang and Yu-Hong Dai. A framework of constraint preserving update schemes for optimization on stiefel manifold. Mathematical Programming, 153(2):535-577 2015.
Michael Jordan, Tianyi Lin, and Emmanouil-Vasiln for min-max optimization in geodesic metric spaces. Advances in Neural Information Processing Systems, 65744, 2022.
KaaakisGGaragkounii. First-order algorithms
Hiroyuki Kasai, Hiroyuki Sato, and Bamdev Mishra. Riemannian stochastic recursive gradient algorithm. In International conference on machine learning, pp. 2516-2526 2524. PMLR, 2018.
Aparajita Khan and Pradipta Maji. Multi-manifold optimization for multi-view subspace clustering. IEEE Transactions on Neural Networks and Learning Systems, 33(8):3895-3977, 2021.
Max Kochurov, Rasul Karimov, and Serge Kozlukov. Geoopt: Riemannian optimization in pytorch.
arXiv preprint arXiv:2005.02819 2020.
John M Lee. Introduction to Riemannian manifolds, volume 2. Springer, 2018.
John M Lee and John M Lee. Submanifolds. Introduction to smooth manifolds, pp. 188, 2012.
Jing Li, Quanxue Gao, Qianqian Wang, Cheng Deng, and Deyan Xie. Label learning method based on tensor projection. In Proceedings of the 30th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, pp. 1599-1609, 2024a.
Jun Li, Li Fuxin, and Sinisa Todorovic. Efficient riemannian optimization on the stiefel manifold via the cayley transform. arXiv preprint 200X::00022.//13 2020.
Shuai Li, Yingjie Zhang. Hongtu Zhu, Christina Wang. Hai Shu, Ziqi Chen, Zhuoran Sun, and Yanfeng Yang. K-nearest-neighhor local sampling based conditional independence testing. Advances in Neural Information Processing Systems, 36, 2024b.
Yeqing Li, Feiping Nie, Heng Huang, and Junzhou Huang. Large-scale multi-view spectral clustering via bipartite graph. In Proceedings of the AAAI conference on artificial intelligence, volume 29, 2015.
Hongfu Liu, Ziming Huang, Qi Chen, Minggin Li, Yun Fu, and Lintao Zhang, Fast clustering with flexible balance constraints, In 2018 IEEE International Conference on Big Data (Big Data), pp. 743-750, 2018.
13

---

Published as a conference paper at ICLR 2026
Alexander Ly, Maarten Marsman, Josine Verhagen, Raoul PPP Grasman, and Eric-Jan Wagenmakers. A tutorial on fisher information. Journal of Mathematical Psychology, 80:40-55. 2017.
Peter Macgregor. Fast and simple spectral clustering in theory and practice. Advances in Neural Information Processing Systems, 36, 2024.
Peter Macgregor and He Sun. A tighter amalysis of spectral clustering, and beyond. In International Conference on Machine Learning, pp. 14717-14742. PMLR, 2022.
Mayank Meghwanshi, Pratik Jawanpuria, Anoop Kunchukuttan, Hiroyuki Kasai, and Bamdev Mishra. Mctorch, a manifold optimization library for deep learning. arXiv preprint 1.1.01181,,0111,, 2018.
Bamdev Mishra, NTV Satyadev, Hiroyuki Kasai, and Pratik Jawanpuria. Manifold optimization for non-linear optimal transport problems. arXiv preprint arXiv:2103.00902, 2021.
Yujie Mo, Zhihe Lu, Runpeng Yu, Xiaofeng Zhu, and Xinchao Wang. Revisiting self-supervised heterogeneous graph learning from spectral clustering perspective. Advances in Neural Information Processing Systems, 17:43133-4311 2025.
Anindya Mondal, Jhony H Giraldo, Thierry Bouwmans, Ananda Chowdhury, al. Moving object detection for event-based vision using graph spectral clustering. In Proceedings of the IEEF/CVF international conference on computer vision, pp. 876-884, 2021.
Andrew Ng, Michael Jordan, and Yair Weiss. On spectral clustering: Analysis and an algorithm. Advances in neural information processing systems, 14, 2001.
Du Nguyen. Closed-form geodesics and optimization for riemannian logarithms of stiefel and flag manifolds. Journal of Optimization Theory and Applications, 194(1))142-166, 2022.
Feiping Nie, Chris Ding, Dijun Luo, and Heng Huang. Improved minmax cut graph clustering with nonnegative relaxation. In Machine Learning and Knowledge Discovery in Databases: European Conference, ECML PKDD 2010, Barcelona, Spain, September 20-24, 2010, Proceedings, Part II 21, pp. 451-466. Springer, 2010.
Feiping Nie, Xiaoqian Wang, and Heng Huang. Clustering and projected clustering with adaptive neighbors. In Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining, pp. 977-986, 2014.
Feiping Nie, Xiaoqian Wang, Michael Jordan, and Heng Huang. The constrained laplacian rank algo- rithm for graph-based clustering. In Proceedings of the AAAI conference on artificial intelligence, volume 30, 2016.
Feiping Nie, Jingjing Xue, Danyang Wu, Rong Wang. Hui Li, and Xuelong Li. Coordinate descent method for & k-means. IEEE Transactions on Pattern Analysis and Machine Intelligence, 44(5): 2371-2385, 2021.
Feiping Nie, Huimin Chen, Heng Huang, Chris HQ Ding, and Xuelong Li. Learning a subspace and clustering simultaneously with manifold regularized nonnegative matrix factorization. Guidance, Navigation and Control, 2024.
Michael Overton and Robert S Womersley. Second derivatives for optimizing eigenvalues of symmetric matrices. SIAM Journal Matrix Analysis and Applications. 16(3):697-718, 1995.
Peter Petersen. Riemannian geometry, volume 171. Springer, 2006.
Yazhou Ren, Jingyu Pu, Zhimeng Yang, Jie Xu, Guofeng Li, Xiaorong Pu, S Yu Philip. and Lifang He. Deep clustering: A comprehensive survey. IEEE transactions on neural networks and learning systems, 2024.
Jorma Rissanen. Fisher information and stochastic complexity. IEEE transactions on information theory, 42(1):40-477 1996.
14

---

Published as a conference paper at ICLR 2026
Meitar Ronen, Shahaf E Finder, and Oren Freifeld. Deepdpm: Deep clustering with unknown number of clusters. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 9861-9870, 2022.
Farid Saberi-Movahed, Kamal Berahman, Razieh Sheikhpour, Yuefeng Li, and Shirui Pan. Nonnega- tive matrix factorization in dimensionality reduction: A survey. arXiv preprini arXiv:2405.03615 2024.
Hiroyuki Sato. Riemannian optimization and its applications, volume 670. Springer, 2021.
Hiroyuki Sato. Riemannian conjugate gradient methods: General framework and specific algorithms with convergence analyses. SIAM Journal on Optimization, 32(4):2690-2711, 2022.
Martin JA Schuetz, J Kyle Brubaker, and Helmut G Katzgraber. Combinatorial optimization with physics-inspieed graph neural networks. Nature Machine Intelligence, 4(4):367-377, 2022.
Han Shen and Tianyi Chen. On penalty-based bilevel gradient descent method. In International Conference on Machine Learning, pp. 30992-31015. PMLR, 2023.
Dai Shi. Junbin Gao, Xia Hong. ST Boris Choy, and Zhiyong Wang. Coupling matrix manifolds assisted optimization for optimal transport problems. Machine Learning, 110:533-558, 2021.
Liangliang Shi, Zhaoqi Shen, and Junchi Yan. nouble-bounder optimal transport for advanced clustering and classification. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 38, pp. 14982-14990, 2024.
Nimita Shinde, Vishnu Narayanan, and James Saunderson. Memory-efficient approximation al- gorithms for max-k-cut and correlation clustering. Advances in Neural Information Processing Systems, 34:8269-888 8281,,28,,2000 2021.
Oleg Smirnov. Tensorflow riemopt: a library for optimization on riemannian manifolds. arXiv preprint arXiv:2105.13921, 2021.
Siti Noraini Sulaiman and Nor Ashidi Mat Isa. Adaptive fuzzy-k-means clustering algorithm for image segmentation. IEEE Transactions on Consumer Electronics, 56(4))2661-2666 2266-268, 2010.
Yanfeng Sun, Junbin Gao, Xia Hong, Bamdev Mishra, and Baocai Yin. Heterogeneous tensor decomposition for clustering via manifold optimization. IEEE transactions on pattern analysis and machine intelligence, 38(3):476-489, 2015.
Yue Sun, Nicolas Flammarion, and Maryam Fazel. Escaping from saddle points on riemannian manifolds. Advances in Neural Information Processing Systems, 32, 2019.
Shingo Takemoto, Kazuki Naganuma, and Shunsuke Ono. Graph spatio-spectral total variation model for hyperspectral image denoising. IEEE Geoscience and Remote Sensing Letters, 19:1-5, 2022.
Ryan J Tibshirani. Dykstra's algorithm, admm, and coordinate descent: Connections, insights, and extensions. Advances Neural Information Processing Systems, 30, 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,201 2016.
Nilesh Tripuraneni, Nicolas Flammarion, Francis Bach, and Michael I Jordan. Averaging stochastic gradient descent on riemannian manifolds. In Conference On Learning Theory, pp. 650-687.
PMLR, 2018.
Anton Tsitsulin, John Palowitch, Bryan Perozzi, and Emmanuel Müller. Graph clustering with graph neural networks. Journal of Machine Learning Research, 24(127):1-21,
2023. 2023.
Diego Valsesia, Giulia Fracastoro, and Enrico Magli. Deep graph-convolutional image denoising. IEEE Transactions on Image Processing, 29:8226-82311 2020.
Nate Veldt. Cut- matching games for generalized hypergraph ratio cuts. In Proceedings of the ACM Web Conference 2023, pp. 694-704, 2023.
15

---

Published as a conference paper at ICLR 2026
Nisheeth K Vishnoi. Geodesic convex optimization: Differentiation on manifolds, geodesics, and convexity. arXiv preprint arXiv: 2018.
Jie Wang, Haiping Lu, Konstantinos N Plataniotis, and Juwei Lu. Gaussian kernel optimization for pattern classification. Pattern recognition, 42(7):1277-2247, 2009.
Yangtao Wang, Xi Shen, Shell Xu Hu, Yuan Yuan, James I. Crowley, and Dominique Vaufreydaz. Self- supervised transformers for unsupervised object discovery using normalized cut. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp.
2022.
Melanie Weber and Suvrit Sra. Riemannian optimization via frank-wolfe methods. Mathematical Programming, 199(1):525 556, 2023.
Zaiwen Wen and Wotao Yin. A feasible method for optimization with orthogonality constraints. Mathematical Programming. 142(():397-434, 2013.
Fangyuan Xie, Jinghui Yuan, Feiping Nie, and Xuelong Li. Dual-bounded nonlinear optimal transport for size constrained min cut clustering. arXiv preprint arXiv:2501.18143 2025.
Jh Yuan, Zhuo Liu, and Feiping Nie. Riemannian fuzzy k-means on product manifolds, 2026. URL https:/ / opwnreview..e.ee////p net /ttp:/ommmmmd-c//a0170
Jinghui Yuan, Hao Chen, Renwei Luo, and Feiping Nie. A margin-maximizing fine-grained ensemble method. arXiv preprint arXiv:2409.12844 2024a.
Jinghui Yuan, Weijin Jiang. Zhe Cao, Fangyuan Xie, Rong Wang, Feiping Nie, and Yuan Yuan. Achieving more with less: densor-ppimizziion-ppeeed ensemble method. arXiv preprint arXiv:2408.02936, 2024b.
Jinghui Yuan, Chusheng Zeng, Fangyuan Xie, Zhe Cao, Mulin Chen, Rong Wang. Feiping Nie, and Yuan Yuan. Doubly stochastic adaptive neighbors clustering via the marcus mapping. arXiv preprint arXiv:2408.0290 2024c.
Jinghui Yuan, Fangyuan Xie, Feiping Nie, and Xuelong Li. Riemannian optimization on relaxed indicator matrix manifold. 2025. URL https: /atrxvivvv.../////////a///aaaavm
20505.
Alp Yurtsever and Suvrit Sra. Ccep is frank-wolfe in disguise. Advances in Neural Information Processing Systems, 35:35352-35364, 2022.
Chao Zhang, Xiaojun Chen, and Shiqian Ma. A riemannian smoothing steepest descent method for non-lipschitz optimization on embedded submanifolds of n. Mathematics of Operations Research, 49(3):1710-17333 2024.
Xiaowei Zhao, Feiping Nie, Rong Wang, and Xuelong Li. Improving projected fuzzy k-means clustering via robust learning. Neurocomputing, 491:34-43, 2022.
Huasong Zhong, Jianlong Wu, Chong Chen, Jianqiang Huang, Minghua Deng, Liqiang Nie, Zhouchen Lin, and Xian-Sheng Hua. Graph contrastive clustering. In Proceedings of the IEEE/CVF international conference on computer vision, pp. 9224-9233, 2021.
Bingxin Zhou, Ruikun Li, Xuebin Zheng, Yu Guang Wang, and Junbin Gao. Graph denoising with framelet regularizers. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2024.
Shixiang Zhu, Liyan Xie, Minghe Zhang, Rui Gao, and Yao Xie. Distributionally robust weighted k-nearest neighbors. Advances in Neural Information Processing Systems, 01-29088-29100 2022.
Wei Zhu, Feiping Nie, and Xuelong Li. Fast spectral clustering with efficient large graph construction. In 2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp.
2492-2496, 2017. doi: 10. 10.1009CCASSPP22117795566
Xiaojing Zhu. A riemannian conjugate gradient method for optimization the stiefel manifold. Computational optimization and Applications, 67:73-110 10, 2017.
16

---

Published as a conference paper at ICLR 2026
CONTENTS
1
Introduction
1
2 Preliminaries
Riemannian Toolbox
3.1
Definition of the Relaxed Indicator Matrix Manifold
3.2 Riemannian Optimization Toolbox for the RIM Manifold
3.3
Comparison Analysis of Time Complexity
RIM Manifold for Graph Cut
Experiments
5.1
Experimental Setups
7
5.1.1
Experiment Setup
5.1.2
Experiment Setup
5.1.3
Experiment Setup
5.1.4
Experiment Setup
5.2
Experimental Results
5.2.1
Result of Experimental
5.2.2
Result of Experimental
8
5.2.3
Result of Experimental
5.2.4
Result of Experimental
6
Limitations
7
Conclusion
10
Acknowledgement
10
9 Statement
11
Appendices
20
A Proof's of Theorems
20
A.I Proof of Lemma
20
A.2 Proof of Lemma
20
A.3 Proof of Lemma
21
A.4 Proof of Lemma
22
A.5 Proof of Theorem
22
A.6 Proof of Theorem
23
A.7 Proof of Theorem 3
25
17

---

Published as a conference paper at ICLR 2026
A.8 Proof of Theorem 4
26
A.9 Proof of Theorem 5
28
A.10 Proof of Theorem
29
A.11 Proof of Theorem 7
31
A.12 Proof of Theorem
33
B Preliminaries
34
B.1
Notations
34
B.2 Introduction to Riemannian Optimization
................................
35
B.3
Introduction to Related Manifolds
36
B.3.1
Single Stochastic Manifold
................................
36
B.3.2
Doubly Stochastic Manifold
37
................................
B.3.3
Stiefel Manifold
38
B.4
Manifold-based Machine Learning Algorithms
................................
39
B.4.1
Algorithms on the Single Stochastic Manifold
39
B.4.2
Algorithms on the Double Stochastic Manifold
B.4.3
Algorithms on the Steifel Manifold
39
B.5 Other Related Work and Background Introduction
40
C
Optimization Algorithms the RIM Manifold
41
C.1 Gradient Descent on the RIM Manifold
41
C,2 Conjugate Gradient Method on the RIM Manifold
41
C.3 Trust Region Method on the RIM Manifold
41
Details of the Experimental Setup
43
D.I
Experiment 2 Setup
43
D.2
Experiment Setup
43
D.3
Experiment 4 Setup
43
D.3.1
Clustering Accuracy (ACC)
44
D.3.2
Normalized Mutual Information (NMI)
D.3.3
Adjusted Rand Index (ARD)
45
D.3.4
Introduction of Real Datasets
45
D.3.5
How Choose and
46
Additional Experimental Results
47
E.1 Results of Experimental
47
E.2 Results of Experimental 2.
47
E.3 Results of Experimental
49
E.4 Results of Experimental 4
49
18

---

Published as a conference paper at ICLR 2026
F
RIM Manifold Equipped with Fisher Metric
51
E1
Dimension and Tangent Space
51
F2
Riemannian Gradient, Riemannian Connection and Riemannian Hessian
51
E3
Retraction Mapping
52
F.4
Which to Use?...
52
Explanation regarding details
53
G.1 Answer to the Question 1:
53
G.2 Answer to the Question 2:
53
G.3 Answer to the Question 3:
53
Reference Code for RIM Manifold Riemannian Toolbox
54
19

---

Published as a conference paper at ICLR 2026
A PROOFS OF THEOREMS
A.1 PROOF OF LEMMA 1
Our relaxed indicator matrix set M = (X| X1, 1,, <,, <u, 0} forms an embedded submanifold of the Euclidean space, with dim M = (n (n-1)). We refer to it as the Relaxed Indicator Matrix Manifold.
Proof. The set M can be viewed the intersection of three sets: M ( < x In < u, X > 0} = 22 nl2 n 923, where 22 = {X | X > 0, X1, = 1n). Q2 =X | >
1). and Ω2 {X X"1, < uJ. Consider the differential of the local defining function for the set
Q1 i.e.
(X +x++V) - In - (XL - 1n)
tV
- D(X1e-1)[VV = lim
lim t-+
V1.
(23)
r-0
Consider the null space of D - [V], given by Ker(D(X1e In) (V) Vl. = 03.
The dimension of this null space is
nc dim (D( (1001,-1)[111)) n
1)c
(24)
In addition, since S2g = | > 1) and 23 = {X X < u]. take U2 as an example. For any directional vector U, there must exist dy 0 such that (X + òoU)"I, In >> I. Thus, both f2 and G2 are open sets. According to Theorem (Petersen, 2006), O forms a manifold, and Q2 and Ωg are open sets. The intersection of open set with manifold remains manifold. Therefore, M = {X X1. = l, I < 1",, < u, X > 0} = Ω 0 22 n 2g is still a manifold, and dim(M) dim(Ker(D(XI.. I.)(VI) (n 1)e.
We refer to .M as the Relaxed Indicator Matrix manifold, abbreviated as the RIM manifold.
A.2 PROOF OF LEMMA 2
By restricting the Buclidean inner product (U - E' U,yV onto the RIM manifold M, the tangent space of A at X is given by TxM {U U1, 0}. For any function H. its Euclidean gradient is Grad H(F). the Riemannian gradient grad H(F) is expressed as following.
grad H(P) = Grad H(F) Grad H(F)1.17
(25)
Proof. According to the definition of tangent space,
TxM Ker(D(XX,, I.)(UT) (U U),
(26)
Let GradH be the gradient of 7 in the Euclidean space. Then, GradH = GradHl + GradH
where GradH represents the component of GradH parallel to M, and GradH represents the component perpendicular to M.
By the definition of the Riemannian gradient,
DH[V]- (GradH, V) (gradH, V)x, V Txx
(27)
Here, (gradH, V)x denotes the inner product equipped on the manifold at X. When (gradH, V)x coincides with the Euclidean inner product, we have
(GradH, V) (GradH4 V) (GradH V) (GradH4), V) = (gradH,, V)x
(28)
for V E M,, since (GradH<1, V) = 0 for V E M. By the Ritz representation theorem, in this case, gradZ is the orthogonal projection of GradH onto the tangent space. The next step solve
the optimization problem:
min UE(U]C1. (E(U)(l,=o)
L
min
IU |GaddH
(29)
The Lagrangian function for the optimization problem is given by: L =11lU |
Taking the gradient with respect to U, we have:
V L U GradH + al? 0
(30)
20

---

Published as a conference paper at ICLR 2026
Solving for U. obtainU Grad7/ a1.". Since U1, O, substituting U gives GradH/ a = GradH/1c ca = 0, which implies a = 1 GradHLe. Therefore, the Riemannian gradient is following.
gradH
GradH GradH GradH1.4e
(31)
A.3 PROOF OF LEMMA 3
For the manifold M, there exists unique connection that is compatible with the inner product and ensures that the Riemannian Hessian mapping is self-adjoint. This connection is given by following. V is the Riemannian connection in Euclidean space.
VvU'= VV-=-VU--VvVia
(32)
Proof. First, we need to prove that the connection compatible with the inner product, which means proving W (U, V) (V U,V) + (U, VwVV We have the following equation
W V) =
= =
Ve [ ;=1==1
D(U,))))))
n n
V(,,) D(V (U. ((V) + (D(U)W))
(U, D(V)[W] - ((V) + (D(U)[W) D(
C + (U. ( ( V).
(33)
Since the standard inner product in Euclidean space is chosen, we have
(U. (.-D(>)Wn)a = -u(U" D(V))W)aa)) (((D(V))))))))))))
(34)
(35)
The last step equals zero because U E M. which implies that U1 = 0. In the Euclidean space, the connection V U is defined as [V)) Furthermore, we have:
VU, (U, - ((V)(W)11)) + (D(U))) D(UO)w]aa
(36)
= V. - + (V - vwll..
(37)
= (U, + (V, VV,
(38)
The second step is to prove that the Hessian map obtained from the connection is self-adjoint. That is, we need to prove [U, V = V.V - V.U, where [U, Vl is the Lie bracket, and [U, Vlf U(V()) (U)). with / being smooth scalar field the manifold M. U and V are tangent vectors of the RIM manifold M. i.e., U.V E TxM. Let U and V be smooth extensions of U and V in the neighborhood of M, satisfying ÜlM U and VlM V. We have [U, DV(Ū- DŪ[V] = Vo V -V Thus, we can prove that:
[U, V] [U,
(39)
-VVoV - VDŪ)IM
(40)
= =Projm =Projm(Vol - (VoV--VVV) -VP0)IM
(41)
-VoV -VoV- - VyU +
(42)
c
-VOV-Vi
(43)
21

---

Published as a conference paper at ICLR 2026
This equality, (VoV Vpūjim = = -Vpo--Ullla holds because [U, VJ is defined in the tangent space of M. Therefore, the expression (VoV VpUUIM and its projection onto the langent space of M must be equal.
A.4 PROOF OF LEMMA 4
For the manifold M equipped with the connection VvU VvUll the Riemannian Hessian mapping satisfies following.
hess H[V] Hess H[V] Hess HIV]1o16
(44)
Proof. The Riemannian Hessian is defined
UU] Vu grad H = Vu Grad H
Grad Hlol
(45)
Using the definition of the Riemannian connection V, have
hess H(U = V grad H = =D ( Grad H Grad H1..1"
[U]
(46)
Grad H(X + (0) Grad H(X)
Grad H(X + (U)1,1
Grad H(X)1.13
lim
lim f0
(47)
= Hess HVV- Hess H[V]1.12
(48)
A.5 PROOF OF THROREM
Let (xV) argmin FemminFeMMlle ||F 1- +V) X M. Then
reg (X (V) (0, + IV 1,, 1)) (())
(49)
where (t),,w((t),p))) (1) are Lagrange multipliers. Moreover, there exists d such that for E (0, -(()) 1,, Inp(()² O. and the Retraction satisfies the following. Where D
(()=
denotes the Levi-Civita derivative.
d
D (x(())) V, dt Rx 1= = 0
Rx(0)= =
(50)
dt
Thus, Rx(IV) a geodesic.
Proof. First, the Lagrangian dual function of the original problem is follows:
CF,,u,,, p) = ; - F + IV -F - 1,,) - - (F" - u) + p FT (51)
Where v, and p are the corresponding Lagrange multipliers, satisfying 0. w 0, 0. Let 8L
BF = 0, then we have the following formular:
dC
X + +1 + ",wp lnp" -VV
(52)
OF
v1 T с That is, F X IV
- 1, + 1mp Since lies on the manifold M and 2 0, the final
result is:
F" = argmin (X + 12 max (0, + -V--VV-l - -1? I,w" + Inp").
(53)
It be proven that F* satisfies the KKT conditions of the original problem. For different the values of the Lagrange multipliers v, w, p vary, and they are functions of t: v(t), ((t). The next step is to prove the three properties of the second-order Retraction, Rx (0) = X, (x(V))) =
V. R'x (IV) 0. First, consider Rx (0) argmin FF XII2..
22

---

Published as a conference paper at ICLR 2026
Since X M, we have Rx (0) X. Additionally, since F*) max (0, X
tV =1." I,w" + Inp")) l-o max (0, X =(0)13 I,w(0)² + Inp())²). We know that
X = max (o, - - - (,0( 1,)()))).
According to the definition, calculate:
(t) (0)
max (0.X tV =1 1,w" + +1mp") -F*(O) (0)
dt (x = lim [-0
lim 20
(54)
Since X E M. we know that Xij > 0, X1, = Ln and / < X1n < u. Furthermore, since V E M, there exists a ô > 0 such that for E (0, 8). we still have (X +V)); > 0, (X +8V)L I,,
and < (X + I < u. This means that for t E (0, 5), we have Rx (EV) argmin
(X + I and since (X +tV E M, it follows that (xV+ = (t) = X + tV. Therefore,
we have:
d Rx(IV)
(1) (0)
X+ IV - X
lim +00
lim 800
(55)
dt
t
t
For de (tVV,, first consider di R'x (tV) exists an interval (0.8) such that (t) X + tV, and within (0,8). without loss of generality, we can assume that v(t)1² (n)(t)" + Inp(t)"
lime-o 1 ( F (t) F*)))) Since there max (0, + vv2 Inw" + 1np")
de
0, and within this interval, X + tV V. Therefore, H lim } (V V) 0. Thus, the Levi-Civita derivative, compatible with the connection, is D
> O. Thus, within (0,
we have d max (0, X +V - =1 I,w" + Inp")
d "
Zt
di ((V)) = O. This concludes the proof.
A.6 PROOF OF THEOREM 2
X X1, = 1,). S2 {X X"1, > >1), and 2 = [X] 1 < w). The primal problem can be solved using the Dykstras (Tibshirani, 2017:
M Ωt U.g Z2 where Ω1 - (X
Boyle & Dykstra, 1986) algorithm by iteratively projecting onto Q2, f22, and Ωg. Specifically:
Proje, (X) = (Xu, + mh) where is determined by Proja, (X)1e = In-
Projo, (X) and Projes (X) are defined similarly. For example,
[X)
if (X')T1 > 44 (1 (1,1-1Xx11, II),,, X if ( - 4ا
Projo, (x))
(56)
where X is the j-th column of X, and 15 is the j-th element of the column vector 1.
Proof. Consider first the orthogonal projection on Q4 which is to solve the optimization problem:
F = arg min || - where 2 = [X X > 0, X. }. The Lagrange function for this problem, incorporating the equality constraint X1, 1, and the inequality constraint X > 0, is:
C(F, n. θ) = _ - -
(F1. (FF 1 1m)- - DonFF
1,.)
(57)
i.j
where 7) E Rn are Lagrange multipliers for the equality constraints, and O 2 are multipliers for the non-negativity constraints.
Since the constraints are separable row-wise, optimize each row F, independently. The row-wise Lagrangian is (, m/. θ,) = 1I x,1 77: (File 1) - E, FiyF.y, Taking the gradient with respect to Fj and setting it to zero:
F X: 1 - θ = 0 =
F X; + 1 +
(58)
By complementary slackness, O, = 0. If Fij > 0. then O2 = 0, implying F= + m. If
Fly = 0, then X,, + 17, + 1, with O; 2 0, hence X,, + 7 < 0. Thus, the optimal solution is:
FrJ max X,y + m, 0) (X n)+
(59)
23

---

Published as a conference paper at ICLR 2026
The multiplier 77 is determined by the equality constraint Fr = = - Σj=1 (X1; + +1,)+ =1
For the projection onto 92, consider the optimization problem: F argmin |IF x1 where l2 = {X | 1,, > lJ. For each column X², solve: minys || x s.t. (F))"1, > 4
If ( > 1j. the constraint is already satisfied: Projes (X') X
If (x I lg. introduce the Lagrangian:
C(F), = F2 121 + ((F))(( X 0
(60)
Taking the gradient of FS, we have the following:
F X Adn =0
p= =x A1n
(61)
Substitute into the binding constraint (Fi)T1n² = =lj:
(x* + Al In =1 + x - (1, - (x))
(62)
Thus, the projection is:
Projoc (x¹)= xxxx
X + (1 (1) - 1,)1
(63)
n
Combining both cases, we have that
(X), [1 ,,,,),,,
(x,,
if (x)) 1 1/2 1x S 1,
Projo, (X²)
(64)
Similarly, for the projection onto Ug we can follow the same procedure and obtain: if ( < Ujj
Proje, (x*) X ( II,,,,, +X, (X*)*,, Uj.
(65)
where X is the j-th column of X, and Uj is the j-th element of the column vector u.
The ultimate goal is to perform an orthogonal projection onto the intersection of three con- vex sets, 02, Ωg. This can be achieved using the von Neumann iterative projection the- orem. However, the von Neumann iterative projection can only guarantee convergence to Ωi n P2 n Qg. but it does not ensure the orthogonal projection, i.e., the solution to the Retraction problem. To address this, we introduce Dykstras's projection algorithm, which performs linear correction to the Neumann projection algorithm at each step, ensur- ing that it achieves the orthogonal projection onto Q n S2 n lg. The algorithm flowchart for Dykstras's 's projection algorithm for the intersection of d convex sets is shown below.
Algorithm 1: Dykstras's 's Algorithm for Projection onto the Intersection of Convex Sets
Input: Closed convex sets Q2 S2 Pa and point y Raxe
Output: Sequence of iterates u (h) converging to the projection onto Ω2 m O Qd
Initialize (0) y. 21 (0)
(0) 2d 0:
while not converged do
(k) U1 1k)
for = 1 to do
u, (k)
(w:-1 1 + 23 (k-1));;
(k)
Proja,
(k) u
(k)
(k-1) 1) (k).
+ 24
end
l: + 4: + 1;
end
10 return 16 (K)
The algorithm iteratively performs Proje, (-), Proja, (-), and Proja, (-), and at each step, a linear correction using ((x) is applied. This ensures the final result the orthogonal projection onto the intersection Ω Ω O Zzz.
24

---

Published as a conference paper at ICLR 2026
A.7
PROOF OF THEOREM 3
Solving the primal problem is equivalent to solving the following dual problem:
1 1 max(0, X + tV - -v12- wnw² + Inp"))l (v, In) (w,, u) + (p.1) (66)
max L C= 220,220 2
where v, w, and are Lagrange multipliers. The partial derivatives of L with respect to v, w. and known, and gradient ascent can be used solving v, w, and p. Finally, Rx (tV) can be obtained using max(0, X + tV - v1 law² + Inp"). The partial derivatives are following.
0L max(0, X + tV v1
w, Th + 1 1np") In
dv
Oc max(0, X + tV =12² wn + 1np 1np"))" y In u
(67)
dw
OC max (0, X + IV u1? - lnw + 1np ) +
0p
Proof. According to the previous theorem, we know that
F = max (0, + tV - =1 - Inw" + Inp"))
(68)
Substituting F" into the Lagrangian function, we obtain
[(v,w, 0) = max (0, X + + v12 - Inw + 1np²)) - X - +V ||
(69)
2
+ -" max (0, X IV - -v1? Inw" + -1pp=) (3) 01 v²1,,
(70)
+ wT max (0, X + tV -v1? InWT + [np²) )²² T 'In - -
(71)
n" max (0, X + IV - -v17 ln + +1mp" T)
T 1n + p71
(72)
Among the Lagrange multipliers v W. P. we have w 2 0 and 2 O.
If (X + +V - v1?- - 1,,T + 1np") < O, then max (0. X + tV v1 - 1nw + 1np²)
0.
which further leads to
C(v,w.p) 1 | 11 + v"1, w u + p=1 2
(73)
At this point, simple differentiation yields:
a Ow
0 L(v,w,p)=l =1
c(),w,p) = -In
[(u,w,p)= uu
(74)
dv
0p
If (X +V - v1 - wnw²" + Inp=))= 2 O, then max (0,X +tV v1 - Inw- + Inp= = X +
VV -- - 1, + +,pp" is worth noting that nT max (0, X + +IV - -v1 - -I, *1 ll E R
T
is a real number, that is,
R max (0. X + vV - v1 wn=" + 1mp") tmmp") )Ic
(75)
tr (vT max (0, X + +V - v1 Inw T + 1np ())))
(76)
max (0, X + tV -v1? - 1, + +1pp") T) [v1e))
(77)
max (0, X + tV -v1. 1,w= 1np=" ) ,ul?))
(78)
At this point, we have
1
L(v,w,,p)
2 |1 + wnw" - 71 - - 1n) - (w, + (p
(79)
+ {X + tV - v1 - I,, + Inp" 11 + In,² - 1np²")
(80)
2 !1 + -1n - (v, 1n) (w, + (p
1
(81)
(x + 1V, iV,v + ln - Inp") 1l111 + lnw 1, 11?
(82)
25

---

Published as a conference paper at ICLR 2026
2 | + Inw 1n0=11 (p) Ln) - (w,, u)
(83)
(0,1) + (x 1V -1 lapp)
-
(84)
At this point, taking derivatives of the Lagrangian with respect to the multipliers vs w, p. we obtain
Oc
(X + tV - -v1 Inw" Inpp )Lp Ins
dv
OC
(X + IV - -u1 - -Inw" Inp² )
(85)
Ow
0L
= -X + +V - v1? -nnw" + 1xp ) + l.
Op
Finally, by consolidating the two cases, we obtain
0C
Qv max(0, X IV u1 Inw 1np I1 las
Anp")l
0C
max(0, X + tV - =1 In 1np / y -t
(86)
Ow
OC 0p
max(0, X + iV - -V Inw + Inp I + I.
After obtaining the gradient, the dual problem can be solved by a simple dual gradient ascent method. It should be noted that the multipliers and p have non-negative constraints, so projection onto the constraints is needed. Specifically, after each gradient ascent step, w and should be projected onto the non-negative constraint. Once v, w, and p are obtained, F can be derived using
max (0, X + IV v1 n, 1np""
(87)
The algorithm flow is as follows:
Algorithm 2: Dual Gradient Projection Ascent Method
Input: Initial values: Po, wb:
Step size K
Constraints: w 2 0, 2 0
Output: Optimized multipliers: v",w", w., p
Initialize vn. wp. p po:
while not converged do
Compute Gradient:;
8C 0c 0c. Bes Op
Update multipliers:; + v + 8c. 0y
w w ac.
p p 0p
Project onto constraints:; max(0,c4 w); f max (0, p);
11
12 end
return Final values U, U.W. P:
A.8 PROOF OF THROREM 4
The Sinkhom--bas based Retraction is defined as
R (tV) = S(X @ exp(tV @ X)) = diag(p")(X (X o exp(tV e X)) diag(g" @ww)
(88)
where p", q", w* are vectors, exp(-) denotes element-wise exponentiation, and diag(-) converts vector into diagonal matrix. The vectors p" :q", w² obtained by iteratively updating the following
26

---

Published as a conference paper at ICLR 2026
equations:
pl++1) =I, 2 ((X @ exp(tV @ x)) (0 o w(())),
q(t-1) max 1 (X exp(tV @ X))" plet2) wilt) tce) (X o exp(tV xx) plett) o glattl) 1o)
(89)
w((+1) min
u 2
This iterative procedure ensures the mapping onto the RIM manifold. The solution R& (IV) = diag(p") (X @ exp(tV (ttV&.X) XX) diag(g" @w") equivalent to solving the dual-bound optimal transport problem (12) with an entropy regularization parameter of I.
P ΣΣ ==
Rx (VV) frgginreM/ ( (oog((X)))(((V)))) +8 8=1
+5
(Ft log( (Ft)) Fs)) (90)
Proof. Introduce Lagrange multipliers v) E R (for equality F1, = 1n), and X, v R u > 0 (for inequalities FT "n, > 1, Fr In < u). The Lagrangian is:
((F. (F. log(X o exp(tV e (x))) + Σ
(Fi (og(F))))
1.
(91)
(Flo (1- 1T + T (FFTFTTT 1n - uu).
Differentiate with respect to Fú, and set to zero, have
[Vi)
X(,) ) log F.j + - X3 Dj
- -log X1y exp
(92)
Simplify using (Xg(Xi, exp(tVis) (x;)) = log Xij +VVij/ Xxz:
Vi
- -Xij
+ log Fij + mh - X33 + Uj 1= =0
(93)
Xzj
Solve for Fiy:
[V X,,
[Vij
Fij = Xij exp
7 + 2ラ v1 Xis exp
-
(94)
X,
Since À and are positive, introduce the following variable substitutions:
p ( = e² 2 In5 e S 1n
(95)
w
Writing the component-win wise wise form into matrix form, we have the following formula.
F diag(p) (X exp(tV C x)) diag(q @ w).
(96)
To construct the iterative format, we first consider the equality constraints. Substitute F into F1, = In;
diag(p) (X @ exp(tV e X)) diag(g @ w)l. = 1n + diag(p) (X @ exp(tV ( 1n
(97)
Further, we can derive the iterative update formula for the row equality constraints.
p = 1n 2 ((X o cxp(tV @ x)) (g w)) = p (k- = 1 @ ((x o exp(tV @ X)) (k) o w())
(98)
Next, considering the column constraint ET 1, > I, substituting F, we obtain:
diag(g w) (X c exp(tV @ x))" diag(p)" 1, > 1 (g ④ w) O ((X @ exp(tV @ xx)T p) >
(99)
By the complementary slackness condition, obtain: X [(g o w) @ ((X @ exp(tV @ x)) p) -11,
(100)
27

---

Published as a conference paper at ICLR 2026
At this point, we discuss the complementary slackness condition.
[(9 o w) o ((x exp(tV @ x))T #b = 91 =1
[( o ((x o exp(tV @ x))" =4 + + (10 ((x o exp(tt o xx)² c w)
(101)
The element- wise iterative update formula is then derived as follows.
2 ((x ④ exp(tV @ x))" po pow). 10),
4 = max
(102)
(k+1) max (1 110 ((X o exp(tV e xx)² p o w) 1o),
(103)
t = max (10 ((X o exp(tV e ex)" 1
(104)
Considering the column constraint FT u, substituting F, we obtain:
diag(q @uee w) (X o exp(tV @ x))² 1 < = (gow) の ((X o exp(tV ox)) p)
(105)
By the complementary slackness condition for upper bounds:
V3 u (q 2 w) @ ((X @ exp(lV e x))")) m) xx)""))], =0
(106)
This leads to two cases:
{ (g w o ((x exp((VV @ x))² 1) # Ujj Vj = 0 - w/ =1
( (oow) @ ((X exp(tV @ x)) T p)]; ujj = =,, (u (uo ((xxexxp(t o x))"" 2) o oo)).
(107)
The element-wise update rule is then:
wy = min (ue ((Ix ((X0exp((t) exp(eV o (x)xx b) 3
(0) D) ((ttt) ( ((4tt))) 110);
(108)
= (k+1) min (u @ ((XX @ X))2 (1tti) o
(109)
-w min
(uo 0 ((xx exp(LV @ x))² p(tt)) ( ((xt2) .lo)
(110)
The final update formula can be obtained as follows.
p(4+1) = 1, 2 ((X @ exp(tV @ X)) (1(+) o w(())),
((1-1)- = max (10 (x @ exp((VV @ x))²² p(a+t) O wll))) lo). w(++1) min u Ø (X exp((V @ x))² p(0t1) o guttt) .10) )
(111)
It is easy verify that the result derived from Sinkhorn is indeed Retraction (Douik & Hassibi, 2019). It can be seen that the F obtained through the Retraction R2xs (IV) minimizes the inner product with log(X @ exp(tV @ X)) under the entropy regularization coefficient of 1. On one hand, this entropy regularization is introduced merely to facilitate computation via the Sinkhorn theorem. On the other hand, the regularization coefficient being I lacks practical significance. Moreover, this Retraction is not second-order Retraction, making its theoretical justification in terms of convergence properties less rigorous compared to the monimminig Retraction. Therefore, the norm-minimizing Retraction is recommended.
A.9 PROOF OF THBOREM 5
Theorem 9. The loss function for the Ratio Cut is given by HL.(F) = tr LF(FFF))))) Then, the Euclidean gradient of the loss function with respect F is:
GradHL,(F) -
1)
(112)
28

---

Published as a conference paper at ICLR 2026
Given the substitutions (F"F)-1 J and PTLF K, the Euclidean Hessian map for the loss function is:
H, -1V-V---((L1---VVV--->---
+
(113)
LF F FWWN+=JJ[(V*+=>>>
(114)
Proof. Let the objective function be (F) (r(FF where J (FF Apply small perturbation 8 to F, yielding the variation:
,H, tr ((6F") LFF + F" L(OF)J F LFJ
(8F)) J)
(115)
Using the cyclic property of the trace and symmetry (L is symmetric, J is symmetric), we simplify to:
8H, = 2tr (8F" (LFJ FF(T
(116)
Thus, the Euclidean gradient is:
(F)
(117)
Apply the direction V to the gradient and compute the directional derivative:
d Hess GradH, (F + (V) 1=
(118)
dt
Expanding the components:
The derivative of LFJ gives LVJ -
The derivative of -FJKJ yields:
- (FNT:++==VVNNNAAANVV":::"------ F+ FソJ
-VJKJ -
(119)
Combining and simplifying:
HessH, = (LVV - VAN LF+F"LVVJ)
(120)
Further, to obtain the Riemannian gradient and Riemannian Hessian mapping, the Euclidean gradient and Euclidean Hessian mapping from the above expressions can be projected onto the RIM manifold. This allows for the optimization of the Ratio Cut loss function on the RIM manifold.
A. PROOF OF THEBREM 6
Theorem 6. For any graph cut problem expressed as H(F) t(((F"L))(FTWF>") where W any symmetric matrix, the Euclidean gradient GradH(() is bounded, and satisfies:
[G)ad 2 ||Llovn W | ||
(121)
a
a²
where
Mmin (W)-
α
(122)
n
and miin (W) is the smallest singular value of the matrix W. This implies that H(F) is Lipschitz continuous.
Proof.
The spectral norm of the matrix F, which is its largest singular value, satisfies:
||FI = max (F)2 5
| <
(123)
therefore, |F & Vo.
29

---

Published as a conference paper at ICLR 2026
Let F1 be the j-th column of the matrix F. Given the constraint 1',, L, the fi-norm of F satisfies || Di=1 Fiz > L. By the Cauchy-Schwee Schwarz inequality, we have:
F
11
1
| 1|1 S 11
F
(124)
Vn
Next, we estimate a lower bound for the smallest singular value of the matrix FTFF. F. For any unit vector 0 E IR, we have:
12
1
Foll
E
(125)
71
n
Therefore, the smallest singular value of the matrix satisfies:
min (F)
(126)
Vn
Since W is a symmetric matrix, its singular values are the absolute values of its eigenvalues, i.e., o,(W) [A,) Using the singular value incquality for matrix products, we have:
Mmin (F (FWF) 2 miin (F)2 Omin (W)-
(127)
Substituting the previously derived min (F) 2 Vn
Omin (W) min; [Ac(W)] we obtain
1 Omin (W) min (W).
min (FFFF)
(128)
Furthermore, the upper bound for the spectral norm of the inverse matrix can be estimated as:
[(F"WF))~
n
(129)
Fmin WF)
mmin (W)
and the a can be presented as
mmnn (W)12
a
(130)
n1
Using the same proof method as A.9, we provide the gradient expression for the general graph cut objective function as:
GradH(
(131)
and with the above technique, estimate its nuclear upper bound.
For (FFFFWF)) Using the sub-multiplicatiiity of the spectral norm (lABl S ||
[ )):
(LF(F"WF)) "Ile S || |Fll |[( ((F= le
(132)
Substituting the known upper bounds:
[FF(F" -
|
F
=
||
FF
m (FT (FTWF) F)
(133)
< |||e
Fl
72
|F @
|||@ vn
(134)
Omin (W)
Next, we consider the second term posed into four parts, namely:
This term be decom-
||W LF - ||F [|
WFF
(135)
For || F Fles. we have the following inequality:
||Fll < ||Wle
||
||]l Vn.
(136)
30

---

Published as a conference paper at ICLR 2026
For ||FT LF @. we have the following inequality:
| LF'le I ll ||
||
|f
||1|
|I & 11
(137)
Combining our estimates with the previous inequality, we obtain:
(FF(""WE)""((""" LF)(F WFF)
WF)
le
(138)
|WFF [FF 'lo
|| LF
||(F=I WF) '1o
(139)
| vn ||اا
W vn. | n (
min (F WF
O (W)
(140)
W L n/22
W @
|
3/2
(141)
min (W)
a²
In summary, we have
ILllov
W l
L 8/2 3/2
[G)adH(F)]ll
2
(142)
a²
where
Min (W)- a
(143)
n
Since
[GradH([F)]pp S minnn, c) [GGaddH(F)))
(144)
it follows that [[GradH([P)]] p is also bounded.
In partieular, for the Ratio Cut, we know that W the identity matrix. Therefore,
| |Lllev
| | (2
3/2
(FradHm(FF)l 2
(145)
a
a
72
Furthermore, since
grad (F) Grad, H(F) Grad, 111)
(146)
it is clear that grad H., (F') is also bounded. An obvious bound is given by
Ilgrad (F)
||raad H., +
1
(IGrad H, (1))l -11122
(147)
which leads to Iigrad H. (F)
LLlovA
12 a²
|
2
2
+
Vnc
a
(148)
2)
2
||
ni
|1/118888/2
||
(2 +
(149)
a
a
where a
1
A.11 PROOF OF THEOREM 7
Theorem For general graph cut problem expressed as H(F) = (((()))((())))""">>>>>> where W is an arbitrary symmetric matrix, the problem is always Lipschitz smooth. Let the corresponding smoothness Lipschitz constant be Q. When applying Riemannian Gradient Descent (RIMRGD) on the RIM manifold with step size if < z. then H(F) converges to critical point at rate of
O (+)). i.e.
min grad H(F())) 12
2 (H(F(()) H(F*))
(150)
0 OSKST
<(T +
where 7 is the total number of iterations, and H(F*) is the global minimum of H(F).
31

---

Published as a conference paper at ICLR 2026
Proof. For general graph cut problem, similar Theorem A.9, the expression of the Euclidean Hessian mapping can be given.
Hes 2[(VJ - LFJJ WVJKJ
(151)
+
(152)
+ WVFJKJsym(V"WF))
(153)
Where (FTWF)-1- J and LF K, and sym(-) denotes the symmetrization operation.
Similar to the previous discussion, can decompose HessH[[] into multiple parts:
[V]le - C -|WVJK//
(154)
+ +IAA//s(((/"W))JJ///} + || J|lW:sssym""WF))
(155)
|IWFFJKJsm((("W)>))
(156)
So the spectral norm of each part is bounded. is not difficult to prove that the spectral norm of [V) is also bounded. Furthermore, it can be shown that the Riemannian Hessian map [Vs also bounded.
I [[Hess [
+ ||HessH//]]ll
1
(157)
Since Theorem A.5 has already proven that obtain geodesics using Dijkstra's algorithm, in the subsequent proofs, we will directly assume the use of geodesics for the retraction process.
Since the Riemannian Hessian map is bounded, let its upper bound be Q. Using the retraction generated by the geodesic, we can expand the function H as follows:
H(R+(V)) H(F) + (F)) V)
2 I
(158)
In the Riemannian Gradient Descent method on the RIM manifold (RIMRGD), by choosing grad H H(F())), and substituting into the upper bound, we obtain:
H(F(++1) S (H(F(()) x xllgrad H(F(()) + Qx²
2 Igrad (F(()))11
(159)
When the step size S Φ. it simplifies to:
S (((()))))
K grad H(F(k))
(160)
This indicates that each iteration, the function value decreases by at least 7 [| Summing the descent over the first & iterations yields:
grad H(F & H(F(")) H(F(k+11) {(((()) H(F"),,
(161)
where H(F*) is the infimum of H(F). Since the right-hand side is bounded, the series Ei |grad H (F converges, and thus
lim Igrad H(F ((((())] = 0.
(162)
From the inequality above, we obtain:
min [g( 2 H(F*)) OSKST
(163)
x(
which implies convergence rate of O (+).
32

---

Published as a conference paper at ICLR 2026
In addition, since the algorithm in Manopt adopts the Wolfe step size, we further provide conver- gence proof of RIMRGD under the Wolfe step-size scheme. Moreover, based on our experiments, usually yields numerical results consistent with those obtained using the Armijo step size.
Condition 1. Equation (157) shows that the Riemannian Hessian hess bounded. Therefore, we have hess(F) I Q. According Lemma 3.5 (Retraction L-smooth) in (Kasai et al., 2018), there exists > 0 such that
[(x+-1) S f(x) + (grad ((x), s) L1 x++1 = Ro(,), 8 8 T.M.
(164)
Condition 2. We adopt the Wolfe step size, i.e.,
[(x + kd) S [(x) + e x (grad (x)). (grad [(x + kd), d) 2 c2 (grad f(x), d),
(165)
where 0 01 < C2 < 1 hyperparameters.
Condition 3. The Ratio Cut loss is clearly lower bounded (according to the real interpretation of Ratio Cut).
Therefore, according to (Sato, 2021), the algorithm converges to a critical point.
A.12 PROOF THEOREM 8
Theorem 8. Suppose we have an objective function H(F) that isu-strongly geodesically convex and Lipschitz smooth on the doubly-stochastic manifold {X]X > 0, XI, 1,, X = r). After relaxing it to the RIM manifold {_ > 0, I,,4 < X < u), let FT be the optimal solution on the RIM manifold and ES be the optimal solution on the original doubly stochastic manifold. Then we have:
L max (Irr- Z11. ||u rl))
[[ F S O
(166)
Proof. First, convexity is necessary. For a nonconvex problem no method can guarantee finding the global optimum, and moreover any locally optimal point may lie arbitrarily far from the global optimum, so no meaningful distance bound can be established.
Without loss of generality, assume 1 < < u. In this case, Fž is also a feasible point on the RIM manifold, and we view F3 perturbation of FT on the RIM manifold. Therefore, the problem reduces to a sensitivity analysis convex optimization under perturbations. By strong convexity, we have:
H(FF) 2 H(F1) + (Grad, H(FF), F2
FT) +
F.1..
(167)
Rearranging terms gives:
- $ - F1 < H(F)) - H(FF). 2
(168)
Define the magnitude of the perturbation A max((Ir 1||,11u rl)). By sensitivity analysis (Boyd & Vandenberghe, 2004), we obtain:
# | Fy H(F)) H(PT) < L dist(F D)
(169)
where D is the perturbation set induced by relaxing the constraint on T = r, and (ist(FF,, D) denotes the maximum distance from Fi to the perturbation set D. The size of D scales proportionally with 2 (max(|Ir- - Z|I, ||u r This completes the proof.
33

---

Published as a conference paper at ICLR 2026
B PRELIMINARIES
B.1 NOTATIONS
Matrices are denoted by uppercase letters, while vectors are denoted by lowercase letters, Let tr(-) the trace of matrix. 1, denotes an n-dimensional column vector of all ones, and Ind represents x then F E Raxe satisfies the property that each row the set of indicator matrices. If F E Ind"xe, contains exactly one element equal
1, while all others are O. The relaxed indicator matrix set is defined as M = {X] XX.. = la. < X 1,,, < X >> 0}, and we proved it can form a manifold M. M represents the tangent space of M at X. (:::) denotes the Euclidean inner product, while )x denotes the inner product on the manifold at X. H represents the objective function, Grad H denotes the Euclidean gradient of H. and grad H. denotes the Riemannian gradient H. Hess H(F) represents the Euclidean Hessian mapping, while hess H(F)) represents the Riemannian Hessian mapping. Rx denotes the Retraction function at X, which generates a curve passing through X, and Rx (tV) represents curve on the manifold obtained via the Retraction function, satisfying d Rx (0)= V. The connection in Euclidean space is denoted as VvU, while the connection on the manifold is denoted as V.U. The differential mapping is represented as [V] Specifically, geodesic (t) is curve the manifold that extremizes the distance between two points. If di (t) = 0, then 7(t) is geodesic. P represents vector transport, which maps the tangent vector V at point X on the manifold to the tangent space Ty M at another point Y.
di
We have compiled all the symbols used in this paper in Table 7, where their specific meanings explained. Additionally, all Riemannian optimization-related symbols used this paper follow standard conventions in the field and also be referenced in relevant textbooks.
Table 7: Notations.
Notation
Description
Inde
The set of indicater matrices All-ones column vectors size n or Laplacian matrix
L
Lower and upper bounds of the column sum the relaxed indicator matrix, both are o-dimcasional column vectors A that forms manifold
M
lnner product defined in Fuclidean space, mapping two Euclidean voctors to a scalar Inner product defined on the tangent space of J. at X Tangent space
Tx.
manifold A at which The objective function to be optimized
space
N
Grad R(P)
Euclidean gradient H i.e. the gradient emberkling space Riemannian gradient of ' at
grad ((F)
VvU
Riemannian connection of the tangent vector field along in Euclidean space
VyU
Riemannian connection of the lungent vector field along V on the mamifold Riemannian Hessian mapping along tangent vector in Euclidean space Riemannian Hessian mapping along tangent vector V on the manifold A curve on the manifold generated at X along the langent vector EV The derivative (tV)
Hess X[V]
bess Rt V]
(x((V)
((V)
a (()
Levi Civita derivative of a (+))
0, where 2r'(0 means (tV) generates geodesic with parameter!
argminn(>
Returs the minimizer of an optimization problem Linear suhmanifolds that require projection
0g 0g.0g Og X,
The i-th row of matrix X j-th column matrix Orthogonal projection the j-th column of matrix X onto the Retums the maximum
X
Projo, (X")
max (n, min(a,
and
Returns the minimum of and Lagrangian function for solving the optimization problem Frobenius norn of matrix
ll
w00), p())
Lagrange multipliers the optimization problem
de
Partil derivatives L respect v(f), w()),
exp(-)
Element wise expomential function
matrix
ding(.)
Converts a vector into a dingonal matrix The differential mapping of H F along V
Sinkhorn function that outputs doubly stochastic matrix
Maps the tangent vector V point X the manifold the tangent space Ty another point Moore Penrose pauudoinverse of matrix
(r(-)
2
Mleeent-wise division
Hadamard product (element-wise multiplication)
34

---

Published as a conference paper at ICLR 2026
B.2 INTRODUCTION TO RIEMANNIAN OPTIMIZATION
Riemannian optimization optimizes functions over Riemannian manifolds, which are smooth mani- folds equipped with metric that defines distance and angles (Meghwanshi al., 2018). It extends classical optimization to non-Euclidean spaces by replacing the Euclidean gradient with the Rieman- nian gradient and so on. Introduced in the 1990s in control theory and signal processing (Edelman et al., 1998; Overton & Womersley, 1995), it has since been widely adopted in machine learning, computer vision, and data science due to ability handle geometric constraints (Carson al., 2017; Khan & Maji, 2021; Boumal, 2023).
The core idea is to respect the manifold's geometry during optimization. Unlike classical methods that assume Euclidean space, Riemannian optimization accounts for curvature. Early methods used steepest descent, while later developments introduced second-order methods like Riemannian conjugate gradient and Newton methods for faster convergence. Recent advancements have expanded this framework to more complex manifolds, such as Stiefel manifold.
The main advantage of Riemannian optimization lies in its ability to perform optimization directly on the manifold, ensuring that the constraints inherent to the problem are naturally respected. For example, in low-rank matrix factorization, the optimization occurs the Stiefel manifold St = {X E Ba*k x = }g}. where I is the identity matrix of size k, naturally respecting the orthogonality constraints of the factor matrices.
In Riemannian submanifold of Euclidean space, the Riemannian gradient grad?((F) F) at a point F E M is defined as the projection of the Euclidean gradient onto the tangent space of the manifold:
gFad = Projorr (Fad)((F))
(170)
This ensures that the optimization process stays within the manifold, preserving its geometric structure.
To solve optimization problems efficiently on manifolds, key operations include the Riemannian gradient, which is used in gradient- based methods. The gradient descent update rule is:
k ( (-a grad))(F(()))
(171)
where Rs is the Retraction map, and 0g is the step size at iteration k. The purpose of the Retraction to update along a curve in the manifold in a specified direction.
For second-order optimization, the Riemannian Hessian hessH(()) is needed. The Hessian captures the curvature of the manifold and provides more information about the local behavior of the function. The Riemannian Hessian is defined
(Fss = V grad(((FF)
(172)
for any tangent vector V E Tr.,, and is used in more sophisticated optimization algorithms to accelerate convergence.
A geodesic is a curve that connects two points on a manifold with an extremal distance, are also important in Riemannian optimization. They are used to guide the optimization process along the manifold and are defined by the differential equation:
d²
dia (t) + [(r(t), x(t) =
(173)
where I are the Christoffel symbols that encode the manifold's curvature (Boumal, 2014; Smirnov, 2021).
The Retraction map Rx (IV) is used to map from the tangent space back onto the manifold after each iteration. A common Retraction map is the exponential map (Kochurov al., 2020; Sun et al., 2019), which can generate a geodesic.
Riemannian optimization (Yuan al., 2026) efficiently handles manifold structures, avoiding artificial constraints and leading to faster algorithms. Second-order methods like Riemannian conjugate gradient (RCG) and Newton methods further improve convergence by utilizing curvature information. The approach versatile, extending to manifolds such as the Stiefel, Grassmannian, and the Relaxed Indicator Matrix (RIM) manifold, which generalizes both single and double stochastic manifolds.
35

---

Published as a conference paper at ICLR 2026
Overall, Riemannian optimization has become crucial tool in solving large-scale, constrained optimization problems, particularly in machine learning, computer vision, and robotics, due to its ability to manage manifold-valuet data and complex constraints.
B.3 INTRODUCTION TO RELATED MANIFOLDS
In this section, we will introduce the single stochastic manifold, the doubly stochastic manifold, and the Stiefel manifold. For each of these manifolds, will provide their basic definitions and discuss optimization methods on these manifolds.
B.3.1 SINGLE STOCHASTIC MANIFOLD
The single stochastic manifold (Sun et al., 2015; Saberi-Movahed al., 2024) consists of matrices where each element is greater than zero and the row sums are equal to one, denoted as (X X O, X1e = ]. with a dimension (n 1)c. The tangent space of a manifold M at a point X is given by Tx.M= IU X1, 0}.
current research, the Fisher information metric typically used as the inner product on the single stochastic manifold M. and is defined as:
Vi Xij
<U. ΣΣ
VU. TxMM,X
(174)
The Riemannian gradient grad (FFF) is the projection of the Euclidean gradient Grad H(F):
grad (F) Projr,.M (Grad (FF) @ F)
(175)
where Projt, M is the projection operator that projects vectors from the Euclidean space onto T.
Specifically, the projection is given by:
Proj (Z) = - (a1)) o a Z1, E
(176)
This projection operation involves matrix multiplication and element- wise operations, with com- plexity of (nc).
In the single stochastic manifold, the Retraction mapping Rx((V) is defined as:
X Rx (tV) = (X @ exp (tV @ X)) @ (x @ exp (V Ø x)),1?),
where the operation ( denotes element-wise multiplication, and @ denotes element-win wise division. The time complexity of this operation involves clement-wise computation and normalization, resulting
complexity (nc).
In the embedded space, the connection is considered with the Fisher metric the set [X]X > 0). According to the Koszul formula theorem, the unique connection the embedded space is given by:
VoV = DV[ 1 ( e
(177)
Based on this, the unique connection the manifold that makes the Riemannian Hessian mapping self adjoint is:
V Projiy. (VoV) - Projiym -
1
(178)
When involving directional derivatives and projections, the complexity of the operation O(nc).
By computing the connection of the Riemannian gradient, one can obtain the Riemannian Hessian mapping on the manifold. The Riemannian Hessian hess H(F)[V] is
1 hess H(F))V] = Projt,M (D grad H(F) 2 (V c grad H(F)) @ F
(179)
where the computation of D grad HVF)[V] involves the Euclidean directional derivative:
grad H(P)(V) DGrad H(F)[V] Grad VVF)VV - -DaaVV]))) (180)
where (Grad (((F) F) The time complexity of this computation involves higher-order
derivatives and projections, leading to a complexity of ( Due to the complexity of the computa-
tion, the coefficient in front of ((n)) is large.
36

---

Published as a conference paper at ICLR 2026
B.3.2 DOUBLY STOCHASTIC MANIFOLD
The double stochastic manifold (Shi et al., 2021; Douik & Hassibi, 2019) refers the set of matrices where each element is greater than 0, the row sums equal 1, and the column sums equal r. Specifically, the manifold defined as:
{X]
X 0, x1. = I,, x 1, =r}
(181)
with dimension (n 1)(c 1). In fact, there are requirements for The more general definition is as follows.
{X X > 0, X1. = In,X" 11=.,,r, =r,"",, 1,,3
(182)
where r is a general vector and the last condition ensures consistency of row and column sums. Generally, simply denote it as (181). The tangent space of the manifold M X is:
TxM = IU X1, 0, x", = 0}
(183)
In current research, the Fisher information metric is also used the inner product on the double stochastic manifold M, defined (U,V)x Σ, UnyYa, VU,V TxM, X M. The
Xzy
Riemannian gradient on the double stochastic manifold given by (n=c):
grad H(F) = - (a17 + In 7 L F) о F. a (I- (7 FT²) ) Ins 7= Grad H(F)
(184)
Here, (F)) represents the Moore- Penorse pseudoinverse of an n x n matrix. Since computing the pseudoinverse requires at least ((n²) operations, this method is impractical for large-scale datasets.
The connection on the double stochastic manifold is defined as an embedded manifold, and in the embedding space, the connection is given by VoV = DV[U] 2 O V) Ø X. Further, the connection on the double stochastic manifold is given by Projy Pr (VvV) = ProjTxM (DV[U)
v)
Projrx,M M denotes the projection into the tangent space of the double stochastic manifold. The projection expression is:
Proje -- (a1% + 1,,3) о OX, a= (1- (z (2x In, B 1-1,--aa
(185)
Indeed, the Riemannian Hessian mapping calculation in the referenced literature involves very complex expressions, including pseudoinverses and other operations with a time complexity of Onn)), making it infeasible for large-scale datasets. In contrast, the proposed RIM manifold in this paper simplifies the calculation significantly, reducing the complexity to O(n).
The Riemannian Hessian is computed follows:
hess H(F)(V] ProjixM (8 (6-1(00))FF)
a é (7 - F-- In
8 1,, FTa
= Grad ((F)) O
7 (a11 1n & F
(I -FFT))
(186)
à = ? (7 FnT)
+ 4
1,
8 1"1, VT a F à
= Hess ((FF)V]] F Grad H(F) о V
4- (a17 + 1,
OF F (al7 + 1, 3T ov
c (FVT + +VF)
The Retraction map uses Sinkhorn to obtain the doubly stochastic matrix. The time complexity of optimization on the doubly stochastic manifold is large, with constant term of O The aboved formulas suitable for the case where n c. However, when n o c, the calculation formula differs slightly, but the time complexity remains the same.
37

---

Published as a conference paper at ICLR 2026
B.3.3 STIEFEL MANIFOLD
The Stiefel manifold (Jiang & Dai, 2015; Li et al., 2020; Zhu, 2017) is the of all matrices whose columns are orthonormal, i.e., St(n, {x E R
X X = 7)-
(187)
can be proven that this set satisfies the requirements for a manifold, and the dimension of this manifold is given by:
dim dim(St(n, c)) = nc </c + 1) At X E St, the tangent space of the Stiefel manifold is given by:
(188)
2
TxSt = 2z | 2²'x + x 0}. Since the Stiefel manifold is embedded submanifold of Raxe its Riemannian inner product is defined as the Euclidean inner product V)x = Eis Uiy Vig.
(189)
The projection operator onto the tangent space Tx is given by:
Proj (W - W
(190)
W -x - 1x(XXZX7))
Based on this, the Riemannian gradient can be directly obtained by projecting the gradient.
grad HFF) Proir H(F)) (W WT),, W Grad G(F))((F)F - F(F" -1a Grad H(FF)"))
(191)
To compute the Retraction the Steifel manifold, the Cayley transform method is used, given by:
w) - (1+) w))
(()) -1-
(192)
Where W = V - WT, is the length on the curve. However, the inversion of (1 - | w) is computationally expensive. To address this, Li et al. (2020) further attempts to use an iterative approach to find the solution. The Retraction is obtained by iteratively solving the following equation:
Y(a) X W(X +Y(a))
(193)
Even so, each iteration still requires multiple matrix multiplications, resulting in a relatively high computational cost.
To obtain the momentum gradient descent on the Riemannian manifold. it is necessary to define the vector transport, which moves tangent vector Vı Tx, St from the Steifel manifold at X1 to the tangent space Tx,St St at X2. This transport operation is denoted as:
P St, - St, VV1 E
E St.
(194)
In fact, this transport operation general in its definition for manifolds. For the Relaxed Indicator Matrix (RIM) manifold, T M Tx,M M for all X1, X2 M, which means that the vector transport is simply P(Vi) = V1 in the RIM manifold. However, this property does not hold on the Steifel manifold. The transport formula on the Steifel manifold is given by:
P(Vi) sr (W - (W-WT)) Wx))X2,
(195)
where V = V2 (X X2), ensuring that the vector is properly projected into the tangent space X2. This projection step ensures the transfer of the vector V1 from the tangent space X, to the tangent space at X2 the Steifel manifold.
As for the computation of the connection and the Riemannian mapping matrix, although the literature does not provide explicit expressions, it can be proven that the expressions for the connection and Hessian map are as follows:
V = Proj se (DV(U))
(196)
hessH((F)[V]] Proitx, (F))))
Using the above Riemannian toolbox, Riemannian optimization can be performed on the Steifel manifold. If the closed-form solution for the Retraction is directly computed. the time complexity On)c). However, by using an iterative approach, the time complexity can be reduced to a large constant factor of On²))
38

---

Published as a conference paper at ICLR 2026
B.4 MANIFOLD-BANE BASED MACHINE LEARNING ALGORITHMS
In this section, we will introduce some classical machine learning algorithms defined on the Single stochastic, Double stochastic, and Steifel manifolds. In general, we assume the data matrix is Z, where Z E R= with n samples and k features. Each row of Z represents sample, and 2 denotes the i-th row of Z.
B.4.1 ALGORITHMM ON THE SINGLE STOCHASTIC MANIFOLD
Fuzzy K-means (Fuzzy C-means, FCM) (Sulaiman & Isa, 2010) is an extension of the traditional K-means algorithm that allows data points to belong to multiple clusters with degrees of membership, rather than being strictly assigned to a single cluster. The core idea is to describe the relationship between data points and clusters through membership matrix, which is suitable for clustering data with fuzzy boundaries.
Let the number of clusters be c, and the membership matrix U E Re where Waj represents the membership degree of the j-th data point in the i-th cluster. The cluster centers are denoted as C {c1, C2,,.., Ce). The optimization goal is to minimize the following objective function:
J(U,
= (U, 2 call
(197)
The constraints are that the sum of the membership degrees for each data point equals I: E May 1 (Vj = 1, 2, ..., r1), and the membership degrees are non-negative: usy E [0, 1]. Where m is the fuzziness coefficient, which controls the degree of fuzziness in the clustering and || call is the Euclidean distance between data point 2j and cluster center c,. Thus, the final objective function and constraints can be written as:
min C(U,C) S.L. UEEXER
0, 01 = E
(198)
This optimization problem is defined over the Cartesian product of the single stochastic manifold and the Euclidean space, which still constitutes a form of a single stochastic manifold.
B.4.2 ALGORITHMS ON THE DOUBLE STOCHASTIC MANIFOLD
ANCMM (Yuan et al., 2024c) is a method for solving constrained problems the double stochastic manifold, which can achieve adaptive neighbor clustering. Its objective function given by:
min SEROXM
Σ
21 - 2,112 S. S||S
(199)
s.t. S 1, S 81 S 1, S
ST, rank(Ls) = n
where S the similarity matrix, and Sij represents the similarity between the i-th and j-th samples. The constraint can be written as:
[X E RR (Ram
X1, = 1, 1
0) n (X E Rxx x Is (200)
where Lg is the Laplacian matrix corresponding to S, and Lg = n c implies that the learned S is naturally c-connected, leading to clusters. Thus, this problem can be viewed constrained optimization problem on the double stochastic manifold.
B.4.3 ALGORITHMS ON THE STEIFEL MANIFOLD
The Min Cut (Fox et al., 2023) is classic clustering method on the Steifel manifold, and its objective function and constraints are given by:
min F LF), F.F.F.F E (F E T F=
(201)
This optimization problem can be solved through eigenvalue decomposition. However, it requires approximately O(n²c) time complexity, and eigenvalue decomposition alone does not provide clustering results. Additional post- processog,,,,gg such as using k-means, is required. Similarly, the derived classic methods such as Ratio Cut and Normalized Cut also classic machine learning
39

---

Published as a conference paper at ICLR 2026
algorithms the Steifel manifold. The expressions for Ratio Cut and Normalized Cut follows:
minf tn FF(F"F)")) st.F E F 11 mine tr L(FFLF(PTDF)))
(202)
S.L. (F Taxe FFFF=
In addition, algorithms such as MinMax Cut (Nie et al., 2010), Principal Component Analysis (PCA) (Abdi & Williams, 2010), Robust PCA (Hubert et al., 2005), and others are also classic machine learning algorithms defined on the Steifel manifold.
B.5 OTHER RELATED WORK AND BACKGROUND INTRODUCTION
In this section, first review our contributions and then introduce other related work beyond
manifold optimization.
As mentioned in our paper, there are currently three main approaches to relaxing the indicator matrix (ours being the fourth). For the first three, the optimization methods themselves have seen little change, but have instead been applied to different models. For example:
The earliest approach relaxes to the singly stochastic manifold (Bezdek et al., 1979), which actually has history of more than 45 years. More recent applications in clustering include (Bao et al., 2024), which employs momentum methods to solve the constraint, and Zhao et al. (2022), which introduces auxiliary variables and updates via coordinate descent. The main drawback of this relaxation is its inability to incorporate prior information about class sizes into the model.
Another line of work relaxes the Stiefel manifold, starting from (Ng et al., 2001), which spurred the development of spectral graph theory and has now history of about 20 years. The basic idea is construct forms like tr tr(1" LF)) and perform spectral decomposition, as in (He al., 2025). The limitation here is that the resulting I'' lacks the interpretability of an indicator matrix, requiring subsequent K-Means step, with a computational complexity of O(n Moreover, this approach also cannot incorporate any class-related information.
more recent direction doubly stochastic relaxation, with representative work Fettal al. (2024), which solves the problem via optimal transport, and Douik & Hassibi (2019), which adopts manifold optimization. The challenge here that the constraints can be overly strict and counterproductive to the model, and manifold optimization still requires Onz)).
Some works in optimal transport are also related to ours. For example, Chapel et al. (2020) introduces Partial Optimal Transport, which is less strict form of optimal transport. This idea is similar to ours in spirit; however, our algorithm designed for arbitrary functions defined on manifolds, whereas theirs focuses classical linear problems.
In addition, Benamou et al. (2014) shows that optimal transport problems can be solved using Bregman Projections. This is close in spirit to the original motivation behind our Retraction design. We further demonstrate that our Retraction corresponds to a geodesic, while also simplifying the overall algorithmic procedure.
40

---

Published as a conference paper at ICLR 2026
C OPTIMIZATION ALGORITHMS ON THE RIM MANIFOLD
In this section, we will introduce three renowned Riemannian optimization algorithms that are utilized in this paper: the Riemannian Gradient Descent method, the Riemannian Conjugate Gradient method, and the Riemannian Trust-Region method. For each algorithm, we will present its fundamental concepts and provide pseudocode. For detailed implementations of these algorithms, one may refer to the sounce manifold optimization package, Manopt (Boumal et al., 2014).
C.I GRADIENT DESCENT THE RIM MANIFOLD
The Gradient Descent on the RIM Manifold method generalizes the classical gradient descent in Euclidean space to Riemannian manifolds by replacing the traditional gradient with the Riemannian gradient, ensuring that the iterations remain on the manifold. The key idea is to utilize the manifold's geometric structure adjust the gradient direction, and then use Retraction to map the updated point back onto the manifold. The process begins with initialization, where an initial point Fo is chosen on the manifold, and step size chosen. In the next step, the Euclidean gradient of the objective function is computed the current point ((k) Then, the Euclidean gradient is projected onto the tangent space of the manifold to obtain the Riemannian gradient, which involves adjusting the gradient by subtracting the normal component. The updated point is then computed along the Riemannian gradient direction, and Retraction (such as exponential mapping or projection) is used to ensure that the new point remains on the manifold. The process continues iteratively until the gradient norm or the change in the objective function becomes smaller than a predefined threshold. The reference pseudo code is in Algorithm
C.2 CONJUGATE GRADIENT METHOD ON THE RIM MANIFOLD
The Conjugate Gradient Method the RIM Manifold introduces conjugate directions to reduce the redundancy in search directions during iterations, thereby speeding up convergence by incorporating information from previous search directions. The core idea is define and update conjugate directions on the manifold. The method begins with initialization, where the initial point Fo is chosen, the initial Riemannian gradient go is computed, and the initial search direction is set do = -9o- Then, the optimal step size in the direction of dy is determined through a line search, using conditions like Armijo's 's rule. The point updated along dg. and Retraction applied to map back onto the manifold. In the next step, the conjugate direction is updated using the current gradient Gk+1 and the previous direction dg, with formulas such as the Polak-Ribière method to compute the new conjugate direction 1 On the RIM manifold, the transport of tangent vectors is equivalent to the vectors themselves. This property simplifies the process the Riemannian Conjugate Gradient Method. The process is repeated until convergence is achieved. The reference pseudo code is in Algorithm 4.
C.3 TRUST REGION METHOD ON THE RIM MANIFOLD
The Trust Region Method on the RIM Manifold constructs local quadratic model in each iteration and constrains the step size within trust region ensure stability. The trust region radius is dynamically adjusted to balance the accuracy of the model with the step size. The method starts with initialization, where the initial point Fo and trust region radius Ag are set. The Riemannian gradient 9k and the approximate Hessian Hş computed at F()) The next step involves solving the constrained quadratic optimization problem in the tangent space, given by:
1 dd Had 2
min dET SI
gh d
(203)
Following this, the method updates the point and adjusts the trust region radius AS based on the ratio of the actual decrease in the objective function to the model's 's predicted decrease. Finally, Retraction is used to project the updated point back onto the manifold. This method known for its strong stability and is particularly suited for highly nonlinear problems. However, it requires frequent Hessian
41

---

Published as a conference paper at ICLR 2026
calculations, resulting in high computational cost. The reference pseudo code is in Algorithm 5.
Algorithm 3: Riemannian Gradient Descent Algorithm RIM Manifold
Input: RIM manifold M (X X1, 1, < X"1, u, 0} Objective function P)), Retraction Rx (tV), transport P. Initial point Fg M
Output: Sequence of iterates { {F()) converging to stationary point of H Initialize k while not converged do
Compute Euclidean gradient Grad H(F())
Compute Riemannian gradient: grad H ((F(k)) Grad H(((())) Grad H(F())111)1
The line search step size: kk(
Perform Retraction: F(k+1)- Ro (1c(0) grad H(F(())
1:
k
end
return F(k)
Algorithm 4: Riemannian Conjugate Gradient Algorithm on RIM Manifold
Input: RIM manifold M = (X X1, I,, < X-1, < u, > 0} Objective function HFF), Retraction Rx (LV), Initial point Fo E M.
Output: Sequence of iterates {F(k) converging stationary point of H.
Initialize k = 0;
Compute initial Riemannian gradient, do grad H(F()));
while not converged do
Compute line search step size ke(k)
Perform Retraction: p(++1) = Pecc (())
Compute new gradient grad H(((k++1)
Compute the conjugate direction d(k+1) grad ((((+)))) (((k)
Compute ((A); g(k) (grad Hy (((p(()))))))))))) ggrad MC ((((()))))))))))))))))))))))))) -grad Hip(hyy
(grad H(F())) H(F()))
l:
k +
10 end
return (())
Algorithm 5: Riemannian Trust Region Algorithm RIM Manifold
Input: RIM manifold M (X XL. Ln, < 1" u, X > 0} Objective function H(F)), Retraction Rx (tV), Initial point Fo E M, Initial trust region radius Ag.
Output: Sequence of iterates [{{)) } converging to stationary point of H
Initialize Initialize Ao while not converged do Compute Riemannian gradient grad H((F(())
Compute the Riemannian Hessian hess ((F(())
Solve the trust region subproblem: A (k) arg minjailaa d H(F(k) + d)
Compute the step size (k) using line search or heuristic method
Perform Retraction: (*++1) R R (1:())))))))())))))
Update the trust region radius A++1;
1c
1: +
end
10 return F())
42

---

Published as a conference paper at ICLR 2026
D DETAILS OF THE EXPERIMENTAL SETUP
D.1 EXPERIMENT 2 SETUP
In the first problem of Experiment 2, due to the particularity the manifold, it is known that the optimal solution manifold M is A, which point the value of the objective function is O. Therefore, by comparing the losses of different algorithms under various parameters, the one with the smallest loss is the optimal result.
D.2 EXPERIMENT 3 SETUP
For Experiment 3, we compared the cases when and o u. For of u, we set / 0.9 and u. = 1.1 LWJ- When = u, the RIM manifold degenerates into the double stochastic manifold, and we can compare it with algorithms on the double stochastic manifold. When = 0.9 [
and 1.1 [ more general methods such as the Trank-Wolfs Algorithm (FWA) and Projected Gradient Descent (PGD) are used for comparison. For the case where the RIM manifold degenerates into the double stochastic manifold, also compared the Riemannian Gradient Descent (DSRGD) and Riemannian Conjugate Gradient (DSRCG) on the double stochastic manifold. brief introduction these algorithms is provided follows:
The Frank- Wolfe algorithm (Xie al., 2025) is a well-known method for solving nonlinear constrained optimization problems. The core idea is to find the direction within the constraint set that is closest the negative gradient direction, and search and descend along this direction to optimize the objective function.
The Projected Gradient Descent algorithm (Chen et al., 2021) is also a method for solving nonlinear constrained problems. The process involves searching along the gradient direction, and when leaving the constraint set, the point is projected back onto the constraint set.
Riemannian optimization on the double stochastic manifold (Douik & Hassibi, 2019): This includes Double Stochastic Riemannian Gradient Descent, Double Stochastic Riemannian Conjugate Gradient methods. The algorithm process is similar to the RIM manifold methods, except that the Retraction and Riemannian gradient computation methods are different.
The PGD method differs greatly from Riemannian optimization methods, including the search direction. Projected Gradient Descent follows the Euclidean gradient, but the Euclidean gradient may contain irrelevant information on the constraint set. Riemannian optimization removes the redundant information and searches along the Riemannian gradient direction.
The Retraction process also differs; the projection process in Projected Gradient Descent may not be easy compute and the result may not be unique, whereas Riemannian optimization choose appropriate Retraction process, which is faster and more convenient.
The generality is also different: Riemannian optimization not only has Riemannian descent but also be naturally extended to methods like Riemannian Conjugate Gradient, Riemannian Coordinate Descent, etc., while Projected Gradient Descent has fewer such extensions.
The convergence properties differ as well; for example, Projected Gradient Descent typically requires convexity to converge the global optimum, while Riemannian optimization only requires geodesic convexity, and there are cases where non-convex problems are geodesically convex.
We compared the final results obtained by optimizing with these algorithms and the total time required, and we organized the data into tables in the main text and appendix, along with visualizations through plotting.
D.3 EXPERIMENT 4 SETUP
For RIMRcut, we apply the same initialization as (Xie et al., 2025) and perform RIM optimization on Reut based on the initialization. When applying the RIM manifold to the Reut, we compare with ten benchmark clustering algorithms across eight real-world datasets. These algorithms include KM-based methods, bipartite graph clustering techniques, and various balanced clustering
43

---

Published as a conference paper at ICLR 2026
approaches. By solving the Ratio Cut problem the RIM manifold, the clustering results more balanced, as the number of samples within each cluster is constrained to a reasonable range. A detailed introduction to each algorithm is provided below.
KM partitions data into predefined clusters by minimizing the sum of squared distances between data points and their corresponding cluster centers. It is simple but sensitive to initial centroids and struggles with non-spherical clusters.
CDKM (Nie et al., 2021) improves KM by utilizing coordinate descent method to directly solve the discrete indicator matrix instead of alternative optimization. It could optimize the solution of KM further.
Reut minimizes the cut between two sets a graph while considering the size of the sets, aiming to balance the partition.
Neut improves on Ratio-Cut by normalizing the cut, balancing the partition while considering the total graph weight. It's better suited for non-convex and unevenly distributed clusters.
Nystrom (Chen et al., 201 1) method approximates large kernel matrices using subset of data, making spectral clustering scalable and efficient for large datasets.
BKNC (Chen al., 2022a) (Balanced K-Means with a Novel Constraint) extends K-Means by introducing a balance-aware regularizer, allowing flexible control over cluster balance. It is solved using an iterative optimization algorithm and achieves better balance and clustering performance than existing balanced K-Means variants.
FCFC (Liu et al., 2018) is an efficient clustering algorithm that combines K-means with balance penalty, ensuring flexible cluster sizes. It scales well large datasets and outperforms existing methods in efficiency and clustering quality.
FSC (Zhu et al., 2017) improves spectral clustering efficiency by using Balanced K-means based Hierarchical K-means (BKHK) to construct anchor-based similarity graph. It achieves high performance on large-scale data.
LSCR (Chen & Cai, 2011) randomly selects landmarks instead of using K-Means, making it faster but potentially less accurate than LSCK in capturing data structure.
LSCK selects representative landmarks via K-Means to construct a smaller graph, reducing computational cost while preserving clustering quality.
To evaluate the clustering performance comprehensively. three metrics are applied, which are clus- tering accuracy (ACC), normalized mutual information (NMI) and adjusted rand index (ARI). The calculation of these three metrics are displayed below.
D.3.1 CLUSTERING ACCURACY (ACC)
Clustering Accuracy (Yuan et al., 2024a;b) measures the proportion of correctly clustered data points by aligning predicted cluster labels with ground truth labels. Since clustering algorithms do not inherently assign specific labels, permutation mapping is applied, often using the Hungarian algorithm, to maximize alignment. The formula for ACC
ACC ((map(N.)
(204)
where 6(a, b) is an indicator function defined as:
8(a, 6) = [1, a 0, otherwise,
(205)
Here, Å is the predicted label, is the true label, n is the total number of data points, and map( is the permutation mapping function that aligns predicted labels with ground truth labels. ACC ranges from 0 to 1, with higher values indicating better clustering performance.
D.3.2 NORMALIZED MUTUAL INFORMATION (NMI)
Normalized Mutual Information (Zhong et al., 2021) quantifies the mutual dependence between clustering results and ground truth labels, normalized to account for differences in label distributions.
44

---

Published as a conference paper at ICLR 2026
It evaluates the overlap between clusters and true classes using information theory. Given predicted partitions {Ci)imi and ground truth partitions {C.]i=1, NMI is calculated as:
n nnn log
Ei=
Eラ= C n
[c, C₁I
NMI
(206)
-== C log [cc
(Ei- |C 1jnnCyyloog -
Here, || |enotes the size of a set, and C₁ CC represents the number of data points belonging to both the i-th predicted cluster and the j-th ground truth class. NMI ranges from Ó to 1, where indicates perfect agreement between clustering results and ground truth. It is particularly effective in scenarios with imbalanced class distributions.
D.3.3 ADJUSTED RAND INDEX (ARI)
The Adjusted Rand Index (Dang et al., 2021) measures the similarity between predicted clustering and ground truth by comparing all pairs of samples and evaluating whether they assigned to the same cluster in both results. A contingency table is first constructed, where each element by represents the number of samples in both predicted cluster C and ground truth cluster Cy. The formula for ARI is:
Eig Σ, + Σ (7) Σ E, (2)
Σ,
(())
ARI(C,C)
2
(207)
where (122)
1اy[8y1111 ARI ranges from 1 to I, where I indicates perfect clustering, 0 represents random assignments, and negative values indicate cossethhnn-randon elustering. ARI is robust to differences in cluster sizes and does not favor large number of clusters.
D.3.4 INTRODUCTION OF REAL DATASETS
The real-world datasets includes: COIL20, Digit. JAFFE, MSRA25, PalmData25, USPS20. Wave- form21 and MnistDatt055. These datasets are selected for their diversity data types (images, waveforms, and biometric data) and their widespread use in benchmarking machine learning and computer vision algorithms, They provide a comprehensive evaluation framework for testing the robustness and generalization capabilities of the proposed methods. The detailed description of them are displayed below.
The COIL20 dataset contains 1,440 images of 20 distinct objects, with each object captured from different angles. Each image has 1,024 dimensions, making it suitable for object recognition and clustering tasks.
The Digit dataset consists of 1,797 instances of handwritten digits, ranging from 0 to 9, Each sample has 64 dimensions, representing low-resolution grayscale images.
The JAFFE dataset includes 213 facial expression images from 10 subjects, covering seven basic emotions. Each image has 1,024 dimensions, making it suitable for facial expression recognition and emotion analysis.
The MSRA25 dataset is a widely used benchmark for face recognition task. It consists of 1,799 grayscale face images, each resized 16x 16 pixels. The dataset includes 12 clusters, representing different individuals or categories.
The PalmData25 2 dataset consists of 2,000 palmprint images, each with 256 dimensions. It includes 100 clusters.
The USPS20 dataset is a subset of the USPS handwritten digit dataset, containing 1,854 instances. Each sample has 256 dimensions, representing grayscale images of digits.
45

---

Published as a conference paper at ICLR 2026
The Waveform21 dataset contains 2,746 instances of synthetic waveform data, each with dimensions. includes clusters.
The MnistDDat005 dataset is subset of the MNIST dataset, containing 3,495 instances of handwritten digits. Each sample has 784 dimensions, representing 28x28 grayscale images. It widely used for digit recognition, classification, and clustering tasks, providing benchmark for evaluating machine learning models.
D.3.5 How TO CHOOSE / AND T
I and u are pivotal parameters within the RIM manifold. When the values of / and are set to be equal, an approximation of the doubly stochastic manifold can be achieved. When I and are not equal, their application to practical problems holds significant meaning, particularly in the context of unbalanced scenarios. For instance, in clustering tasks, the RIM manifold encompasses all indicator matrices, with I and representing the minimum and maximum number of samples within each cluster, respectively. The magnitude of these parameters can be estimated based on the total number of samples and the known number of clusters. Alternatively, they may be assigned according to certain prior knowledge. However, it is noteworthy that in the absence prior information, the values of and be set within broader range. In addition, suitable choice of I and can also be determined through multiple trials.
The parameter in RIM optimization is listed in Table 8.
Table 8: Values of I and on different data sets for RIMRcut
Datascts
COIL20
[0.6*n/c]
[1.2%
Digit
[0.4"n/]]
[1.6%
JAFFB
[0.4%a/c]
[1.6%a/c]
MSRA25
[0.4"u/c]
[1.6%//c]
PalmData25
[0.4*n/c]
[1.8'/e]
USPS20
[0.6%ule]
[2.0%ale]
Waveform21
[0.4*nic]
[188t/e]
MnistData0S
[0.8*n/e]
[1.4%//c]
Subsequently, we will perform clustering using the data in this table and visualize the clustering results, as shown in Figure and Figure 8.
Moreover, we acknowledge that precisely choosing and u is challenging task, as it essentially equivalent to obtaining prior information about the dataset. Our study is conducted under the assumption that such prior information is available. Nevertheless, we also provide a possible way to estimate this prior knowledge, namely by running K-Means" to approximate the cluster proportions. For instance, on the MnistData05 dataset, the estimation yields (1, u) = (0.86 2, 1.22 x "), which is close to the values selected. Among them, it is preferable to select balanced K-Means variant (Chang et al., 2014) or a K-Means method equipped with balance regularization.
COIL20
Digi (0.EX, 1865) (0.4% 16x)
JATTE
MSRA25
FalmDond5 (U.S,, L85x) (ONN,, 1.71%)
USPS20
Wavefomm11
MasDDaaals
(. K-Means
0..1%, 1477)
(022%, L.BX) (0.4%, 16x) (0.4%, 16x)
(0A7S, 132X
00.72x 1.66]]
(086 225)))
(. Selected
(2x6x.12x)
(0.4%, (0.4%,18k)
(0.6%, (0.6x220x)
(04%% 1.4%1188 1.8x)
00.8%, 1440)
Table 9: (I, 2) values from datasets.
At the same time, although the algorithm is sensitive to (1, u), the sensitivity is not high. Taking the MnistData05 dataset as an example, the performance metrics under different values of and u are as follows. Here, ax denotes x
Table 10: Performance under different (1, u) values on the MnistData05 dataset.
[xx,2x)
(0.xx.17x]
(0.5x, 1,5x)
(0.6x,
141x
(0.7x,
1.3x)
(0.8x, 1.1x)
(0.9x, 1.1x]
ACC
61.23
61.61
62.36
63.12
64.26
65.55
66.09
NMI
54.96
55.53
56.68
57.54
58.68
59.35
61.93
ARI
46.02
46.25
49.37
50.73
51.82
52.87
$3,02
In addition, we also provide results showing how the accuracy varies with (1, u) across other datasets in Table 11.
46

---

Published as a conference paper at ICLR 2026
Table 11: ACC metric under different (1, u) values on the other dataset.
(L.M)
[0x,2X] 2x]
(0.3x; 1.7x)
(0.5x,1.6x)
[0.0x, 1.4x)
10.7x, 1.3x)
[0.8x,1.45 1.4x)
0.9x,1.12X
Waseform21
69.16
68.90
69.99
72.17
71.81
74,80
71.08
USPS20
70.39
67.53
68.23
PalmDat225
87.10
$6.15
87.70
86,60
91.05
89.45
91.10
MSRA25
55.86
JAFFE
86.38
92.95
96.71
96.71
96.71
96.71
92.95
Digit
78.46
81,52
84,14
82.47
79.74
80.80
77.90
COIL20
72.92
76.94
77.99
78.13
78.68
77.22
70.49
E ADDITIONAL EXPERIMENTAL RESULTS
In this section, we present additional experimental results to better demonstrate the advantages the proposed algorithm.
E.1 RESULTS OF EXPERIMENTAL 1
Experiment compares the running times of three Retraction methods under different matrix sizes. In Table 12, compare the running times when is not equal The results that run the fastest under each set of experiments are highlighted in red. Additionally, for better visualization, we present three-dimensional bar chart showing the performance of multiple Retraction methods, as illustrated in Figure 3. The experimental results reveal that when the matrix dimension is small, Sinkhorn outperforms the other two methods terms of speed, while Dykstras shows an advantage when the matrix dimension is larger. This conclusion holds true both when equals and when ( does not equal u. While the efficiency of dual method is always inferior than other methods.
Table 12: Table of Execution Time when + for Different Retraction Algorithms(s)
Row&&Cll
1000
Dual
5000
2000
10000
1000
Sinkhom
5000
7000
500
1000
Dylestma
5000
10000
0.004
0.003
0.081
0.044
1.934
2.78
0.005
0.109
0.006
0.010
5.675
0.005
0.149
0.009
0.012
0.100
23.775
0.014
0.018
0.172
15.483
37,027
58,107
03889
0.071
Tinclol
Toreco
100 586
2000
(a) Dual (1 =
(b) Sinkhom (1 w)
(c) Dykstras (1 = u)
Tinetol
100
2000
(d) Dual + Figure 3: Comparison of running time for different Retraction algorithms.
(c) Sinkhorn (1 * u)
(f) Dykstras +
E.2 RESULTS OF EXPERIMENTAL 2
For the first question in Experiment we compare the application of gradient descent, conjugate gradient, and trust-region methods the RIM manifold. The value of cost function and running time of gradient descent and conjugate gradient on RIM manifold are display in Table 13 and 14. As can be from the two tables, regardless of the optimization method employed, the loss function values and running time of the RIM manifold approach are superior to those the doubly stochastic manifold method. This advantage is attributed the lower computational complexity of gradient
47

---

Published as a conference paper at ICLR 2026
Table 13: Cost and Time on the RIM Manifold and Doubly Stochastic Manifold(GGD),
RIM Manifokl
Doubly Stochastic Manifold
Row&Col
Cost
Time
Cost
Time
Size
5000
7000
10000
10000
7030
5
4.74E-14
L114-13
1.05E-13
1.233
0.974
1.225
4.96E-07
6,08E-07
9.01E-07
17.19
18.07
38.73
1.28E-15
4.48E-05
7.04E-15
1.22E-06
7.73E-07
2.39E-06
12.76
19.20
22.45
20
5.39E-14
1.09B-14
1889-11
0.779
1.266
1,914
3.07E-06
2.79E-06
5.46E-06
18.34
20.08
27.02
50
1,95E-15
8.128-14
1.84--13
3.71E-06
6.38E-06
9.27E-06
48.72
37.79
75.39
70
1.72E-13
4.47B-13
1.73E-13
2.350
2.811
4.356
7.91E-06
9.68E-06
1.82E-05
39.37
64.68
56.13
1,58E-15
1.128-14
2.32E-13
1.37E-05
2.99E-05
46.06
93.26
105.8
and Hessian matrix calculations on the RIM manifold. For example, when the matrix size is 100 by 10,000, for the RTR method, the running time is increased by approximately 200 times. For the RGD method, the time required is only one-twenty five of that for the doubly stochastic manifold. As for RCG method, the running time is increased by approximately 75 times. Meanwhile, optimization methods on RIM manifolds often yield solutions closer to zero (the ratio of losses can even reach TE10) compared to methods doubly stochastic manifolds.
The second issue pertains to the problem of image restoration. We introduced varying levels of noise into two images and then compared the visual outcomes of the RIM manifold-based method with those of the DSM-based method in restoring the original images from their noisy counterparts. The visual results are displayed in Figure 4, which also annotates the values the parameter &. Regardless of the intensity of the noise, the images restored by the RIM method are clearer and retain better texture information.
Table 14: Cost and Time on the RIM Manifold and Doubly Stochastic Manifold(RCG).
RIM Manifold
Doubly Stechastic Manifold
Row&Col
Cost
Time
Cost
Time
Size
5000
7000
5
3.74F-14 14
4.63E-13
1.20B-13
0.285
0.375
0.683
8.16F-10
7.78E-10
2.57F-09
4.624
10.22
15.42
10
6,22E-14
2.56E-13
4.92--14
0.161
0.307
0.579
1.81-099
2.71E-09
1,5EE009
11.29
13.96
20
1.03E-13
1.08E-13
1.52E-15
0.396
0.817
0.558
4.75E-09
3.53F-09
2.67E-09
11.91
16.91
17.13
50
5.69E-14
1.56E-13
2.51B-13
0.859
1.047
1.774
3.74E-09
5.49E-09
4.81E009
45.33
70
2.22F-13 13
1.74E-13
6.37B-17
0.932
1.603
1.024
4.36F-09
2.52F-09
4.81E-09
46.83
73.80
60.49
4.03E-09
5.01E.09
7.99E.09
(a) noise=0.3
(b) DSM.{=0.3
(c) RIM.Ç=0.3
(d) noise=0.5
(e) DSM.ç=0.5
(f RIM,=0055
gg) noise=0.9
(h) DSM.
(i) R.M.c=0.7
(i) noise=0.3
(k) DSM.E=0.3
(I) RIM,E=0.3
(m) no5
(n) DSM.
(o) RIM,6=0.3
(p) noisc=0.9
(q) DSM.6=0.7
RIM.E=0.7
Figure 4: mage Denoising Results.
48

---

Published as a conference paper at ICLR 2026
Table 15: Time and Loss Different Optimization Algorithms on Ratio Cut when / o
FWA
PGD
RIMRGD
RIMRCG
RIMRTR
Mettooss&Methoss
Time
Cost
Time
Cost
Time
Cost
Time
Cost
Time
Cost
COIL20
6,220
3,908
6.061
0.7108
8,040
0.5306
2.601
0.494
15.97
0.588
Digit
5.878
0.389
6.063
0.817
7.355
0.652
1.443
0.755
13.92
0.661
JAFFE
0.257
0.207
1.019
0.294
0.116
1.110
0.260
0.154
3.741
0.103
MSRA25
6.238
0.253
6.444
0.048
9.123
0.037
9.787
0.000
15.95
0.033
PalmData25
77.69
16.40
71.73
3,299
25.54
0.984
5.635
6.686
19.05
12.78
USPS20
6.133
1.631
6.109
1,563
7.025
1.544
1.309
1.729
17.72
1.551
Waveform21
12.62
0.405
8.529
0.452
9.650
0.366
1.571
0.373
46.20
0.366
MnistData05
19.86
0.538
17.09
12.36
16.52
1.693
1.876
2.467
30.47
1.677
E.3 RESULTS OF EXPERIMENTAL.
Experiment 3 compared the objective function values and running times of the rim manifold- -aased approach with other solution methods on real datasets when the objective function was the Ratio Cut. The results for the case where I equals u are shown in Table 5, while the results for not equal to u are reported in Table 15. It can be observed that the RIMRCG method achieves the lowest running time on most datasets. Meanwhile, the RIMRGD method can reach the minimum in terms of loss. Furthermore, for each dataset, we have plotted the iteration curves of the objective function values against the number of iterations for various optimization methods. These are displayed in Figures 5 and 6, respectively. From the convergence curves in Figure 5, it is evident that the Rim manifold-based methods enable the objective to decrease more rapidly within a shorter number of iterations. In contrast, the descent curves of the PGD and DSRGD methods more gradual. A similar experimental outcome is also presented in Figure 6.
(a) COIL20
(b) Digit
(c) MSRA25
(d) PalmData25
Figure Comparison of Loss Decrease for Optimization Algorithms on Real Datasets u))
COIL20
Digit
JAFFE
(d) MSRA25
(e) Pam1Data25
(f) USPS20
(g) Waveform21
(h) MnistData05
Figure 6: Comparison of Loss Decrease for Optimization Algorithms Real Datasets (1 + u).
E.4 RESULTS OF EXPERIMENTAL 4
In this section, we mainly provide two supplementary materials. First, we verify whether Riemannian optimization the RIM manifold ensures that the distribution of each column lies within the prescribed range. Second, we visualize the learned indicator matrix to examine whether each entry Fiy E [0, 1] is satisfied.
Figure 7 illustrates the column sum distributions of the relaxed indicator matrix obtained via Rie- mannian gradient descent (RIMRGD) the RIM manifold for different datasets. The dashed lines represent the values of the upper bound u and lower bound As shown, all column sums eventually lie within the specified interval u].
49

---

Published as a conference paper at ICLR 2026
(a) COIL20
(b) Digit
(c) JAFFE
(d) MSRA25
(e) Pam1Data25 Figure 7: Relaxed Indicator Matrix Column Sum Distribution Graph.
(1) USPS20
(g) Waveform21
(h) MnistData05
It is worth noting that, under the specified bounds [1. u/. not all bounds are necessarily active for every dataset. For instance, in the Digit dataset, the lower bound is active, as the sum of the 8th column reaches the lower bound, while column reaches the upper bound u. In contrast, for the MSRA25 dataset, the upper bound is active, but the lower bound is not. For some datasets like COIL20, neither the lower the upper bounds active, possibly because the dataset naturally leads to balanced partitioning.
Pufiallyy
Sonad Lnteeaaee
(a) COIL.20 Gonor
(b) Digit
(c) JAFFE
(d) MSRA25
Chals
(e) PamlData25
USPS20
(g) Waveform21
(h) MnistDataOS
Figure 8: Visualization of the Relaxed Indicator Matrix.
Figure 8 presents the visualization results of the relaxed indicator matrix. It can be observed that each element F.j lies within the interval [o, 1), and the indicator matrix exhibits a clear clustered structure.
This structure indicates distinet clustering results, suggesting that learning the relaxed indicator matrix manifold effectively captures the underlying structure of the graph.
50

---

Published as a conference paper at ICLR 2026
F RIM MANIFOLD EQUIPPED WITH FISHER METRIC
In this section, we will explain why we assign the Euclidean inner product to M = [X]XI, 1n < 1" < u, X > 0} instead of the currently more commonly used Fisher information metric. The RIM manifold is defined as .M- [X| In.& < 1 < u, X > 0}. where the row sums are equal to I. Therefore, each element's 's rows on the RIM manifold can be considered as a probability distribution and can equip Fisher information metric. Specifically, for each point X on the RIM manifold, there exists a tangent space M which is linear space. Previously, the Euclidean metric equipped this linear space, that is, VU, V E TxM, U, >x= D== Eÿ=1 Uny Vij. This section will discuss the impact optimization over the manifold when the Fisher information metric equipped on Tx.M M, that is,
n
Vi
VU.VETxM, U.V
(208)
Xij
To distinguish it from the previous RIM manifold, we call the RIM manifold equipped with the Fisher information metric the Fisher RIM manifold, abbreviated as FRIM manifold.
F.1 DIMENSION AND TANGENT SPACE
Regarding dimension and tangent space, their definitions depend only on the manifold itself and are independent of the metric equipped on it. Therefore, for the set M, whether it is equipped with the Euclidean metric or the Fisher information metric, has the same dimension and tangent space. That is, both the RIM manifold and the FRIM manifold have dimension of (n2 1)c, and the langent space is M = {U UU1 = 0}. The proof can be found in Theorem A.
F.2 RIEMANNIAN GRADIENT, RIEMANNIAN CONNECTION AND RIEMANNIAN HESSIAN
When the Fisher metric is assigned to (X E Rexe > 0}. the gradient of H at X in (X E Raxe| X > 0} is given by GradH
where GradH the Euclidean gradient. At this time, the FRIM manifold is Riemannian embedded submanifold of {X | > 0}. The Riemannian gradient on the FRIM manifold the orthogonal projection under the Fisher metric. The expression of this orthogonal projection is
Proj M (Z) (a1 O a= a= E ER
(209)
The Riemannian connection on the FRIM manifold the orthogonal projection of the connection under the Fisher metric, where the connection on {X E 1 > 0} can be expressed as
VcV= = DV(U) 1 (UO U. E naxe
(210)
The Riemannian connection on the FRIM manifold is given by
V = ProjTxM (V&V) = ProjTxMM (DV(U) - V E Tx E (X E R x 0}
(211)
Furthermore, the Riemannian Hessian mapping is given by
hess H(X)[V] Proji D (grrad H(X)[V] 2 (V grad (odd((X))))))))))) ex)
(212)
can be seen that the Riemannian gradient, Riemannian connection and Riemannian Hessian on the FRIM manifold are the same as those the single stochastic manifold equipped with the Fisher information metric. Since the FRIM manifold itself can be regarded as a Riemannian embedded submanifold of the single stochastic manifold, it is not surprising that these three Riemannian tools are the same as those on the single stochastic manifold. However, the existing Retraction mapping on the single stochastic manifold cannot be applied to the FRIM manifold, because it cannot be guaranteed that the curve generated by the Retraction mapping on the single stochastic manifold will always lie the FRIM manifold (it may violate the column constraint).
51

---

Published as a conference paper at ICLR 2026
F.3 RETRACTION MAPPING
Although the Retraction mapping on the single stochastic manifold cannot be used as the Retraction mapping the FRIM manifold, the Retraction mapping the RIM manifold proposed in this paper can naturally serve as the Retraction mapping on the FRIM manifold. That is, the FRIM manifold naturally has three Retraction methods respectively given by Theorems A.6, A.7 and A.8. However, Theorem A.5 indicates that the result obtained by Theorem A.6 is a geodesic on the RIM manifold. However, on the FRIM manifold, Theorem A.6 is not an orthogonal projection under the Fisher information metric, so the geodesic on the FRIM manifold cannot be obtained. That is to say. although the three methods in Theorems A.6, A.7 and A.8 can all be used as Retractions, none of them is second-order Retraction.
F.4 WHICH TO USE?
Although there is also set of Riemannian tools available on the FRIM manifold, according to the analysis above, the Riemannian toolbox under the RIM manifold and the Riemannian toolbox on the FRIM manifold have almost the same time complexity and can use the same Retraction. However, when using the Dykstras Retraction, a geodesic can be quickly obtained on the RIM manifold, while it is impossible to obtain a geodesic on the FRIM manifold, meaning a second-order Retraction cannot be achieved. This may have a certain impact on the convergence of the algorithm. Therefore, we recommend using the RIM manifold, which restricts the Fuclidean inner product to the manifold rather than the Fisher information metric.
52

---

Published as a conference paper at ICLR 2026
G EXPLANATION REGARDING DETAILS
In this section, we provide detailed answers to three questions raised by the community. These questions essentially represent the special techniques that enable the RIM manifold to operate at singularities.
Question Why is projection chosen as the Retraction onto {X] | X 2 0, X1, = =I, 5
X²1, u What if strict inequalities are required?
Question 2: For the RIM manifold {X X > 0, XI, = L. < 1 < u). does it
become an empty set when = u = r? Why can the algorithm function normally? Question How does behave near the boundary of the RIM manifold?
G.1 ANSWER TO THE QUESTION 1:
Choosing projection as the Retraction is partly due its speed, but more importantly, due to its stability. When strict inequalities are required, the following analysis can be made:
When * u, one can always choose sufficiently small such that the constraints effectively map.
Specifically, given I and u, if we must require the resulting X to satisfy {X X > 0, X1, = 1,, < 1,,, In < u}. we only need to pick an and project onto the set:
{X X 200++ X1, = I,, + x 1n u c}
(213)
Here, can be chosen to be very small, for example, & 10-12 even smaller than machine precision. This is also very easy to implement in the RIM toolbox: simply set u' u and =1=+,, E, and then input (u', (') into the RIM toolbox. At the same time, the exact same method can be used to handle {X X > 0} and {X X > 0 + s}. In fact, when e is sufficiently small, it does not affect the final result at all.
When / = the advantage of choosing projection the Retraction becomes even more apparent, as we will point out in Answer 2.
Furthermore, it must be noted that in scenarios where the RIM manifold is typically used, such as machine learning, one usually does not need to care about whether values are taken at the boundary. In scenarios like classification and clustering, one typically only concerns themselves with the maximum or minimum elements of the Relaxed Indicator Matrix. From a numerical optimization perspective, a correction using can yield significant acceleration, making a worthwhile trade-off.
G.2 ANSWER TO THE QUESTION
When 1 = since the Retraction we actually adopt projects onto {X X 0,X1, =1,,7
X In < uJ. which in this case means projecting onto:
(X X 2, = n, 111,=1)
(214)
This is the reason why the RIM manifold is always able to obtain points on the doubly stochastic manifold when setting =
Moreover, considering the Riemannian gradient, this is equivalent to projecting the Euclidean gradient, causing components that deviate from the tangent direction to become zero. This is equivalent to gradient projection method corrected by a Lagrange multiplier, which accelerates convergence.
G.3 ANSWER TO THE QUESTION 3:
When + u. we can always leverage the equivalence condition mentioned earlier. When considering the properties of the RIM manifold {X X > 0, X1, 1,,/ < 1, < 4} near and u, it is always possible {X X > 0, X1, 1,, I & < XT1,, the boundary, we can always locate extreme points on the boundary, for instance, where X I, = u.
slightly expand the boundaries, i.e., consider the interior properties of + E}. In particular, since the Retraction projects onto
53

---

Published as a conference paper at ICLR 2026
H REFERENCE CODE FOR RIM MANIFOLD RIEMANNIAN TOOLBOX
Function M RIMfactory (n, e, rew, upper, lower)
maxDsiters min (1000, n=c); if sizc (row,
error ('row should be a column vector of size n. ")
end
if size (upper, 1)
error ('upper should be column vector of size c. )
enc
it size (Lower, 1) -=
error ('owerr should be column vector of size c.
end
naname e() sprintf Sdxåd matrices with positive entries Fl_e-i_n,,
<Fi_n<u', n, c);
M. M.dim = ( (n-1)+a;
Mahash (X) ['z' (X()))
M.lincomb @matrixlincomb/
M. zezovee 0 (X) zeros (n, c);
M. M.transp = (x1, x2, d) ProjIcTangent (d)
M. @
U(:);
(x, u) reshape (u, n, c);
M. vecmatarei sometarreimeeeriee ( true;
M. iniener @ @iproduct;
function ip iproduct (X, eta, zeta)
ip sun (eta +z)) end n.norm 0 (x, eta) sqrt (M.inner (8, eta, eta)); M. M.typicaldist () ntc;
rand @random;
function X random (X)
2 - abs (randn (n, c));
X = Dykstras (2, row, lower, upper, maxDSiters);
ena
M. randvec @randomvec;
function eta randomvec (X)
2 randn (n, e) eta ProjToTangent (8);
end M.proj @projection;
function etaproj = projection (X, eta)
etaproj = FronTooangent (cta); end
M. tangent M. proj;
t.tangent2ambient a (X, eta) eta;
M.egrad2rgrad tegrad2rgrad; function rgrad agrad2rgrad (X, egrad)
rgrad ProjToTangent (egrad); enc M.retr @Retraction; function Y Retraction (X, eta, t)
if nargin =
end Y-Dykstras (Xtt row, lower, upper, maxDSiters);
end M.ehess2thess @ehess2rhess;
function rhess chess2rhess (X, egrad, chess, eta) rhess Tooagenng (ehess); end
end
54

---

Published as a conference paper at ICLR 2026
In this section, will provide reference code for the RIM manifold toolbox. Our code is compatible with the well-known open- manifold optimization toolbox Manopt (Boumal et al., 201 14), allowing the direct use of Manopt's algorithms to implement Riemannian optimization on the RIM manifold. The first code block creates a factory named "RIM", which allows for the direct call to the RIM factory obtain the basic description of the RIM manifold, covering the essential information about the manifold and the invocation of basic Riemannian operations.
Dykstras algorithm is one of the methods for implementing Retraction. Its process involves iterative projections and the condition for determining when to exit the loop.
function on [P] Dykstras (M, a, b_1, b_u, if 1_===b_u
tol=le-2;
else
tol-1c-1;
end
rng (1);
[mn, mc] size (M)
p Mj
21 zeros (mn, mc)
z2 zeros (mn, mc);
23 zeros (mn, mc);
for iter - 1:N
for i 1:mn prev_row P(i, :) zi, ;);
P (i, 1) EProjSimplex new (prev_rreee (erev_roow, row, a
z1 (i, :) = prev_row P(i, :);
end
for j 1:me prev_col = P(;, j) + 22(;, ;);
current_sum sum = sum (prev_col); if current_sni >- b_)
22(;, j] 0;
P(; ji prev_col;
else
delta (b (j) surrent_suml / ma;
new_col prev_ccl + delta ones (mn, 1);
22(:, j) prev_ccl new_col;
P(;, j) = new_col;
end
end
for
j 1:me prev_col = P(;, j) 23(;, d); current_ssu = sum (prev_col); if current_sum <= b_u(j)
z3(;, j) 0;
( j) prev_ccl;
e else
delta (o (ÿ) - current_sume / ma;
new col prev_col delta ones (ma, 1);
z3(:, j) prev_col new_col;
P(;, }) = new_col;
end
end
norm (P*ones (ma, 1) -a, 'fro') < tol s8 all >=)))to))) disp ([ (CCovverged at iteration: ", num2str (iter) ]); break;
end
end
end
55

---

Published as a conference paper at ICLR 2026
In the Dykstras algorithm process, the first step is to project onto the simplex, where the projection function is EProjSimplex, _nw. The code for this is provided below. During usage, you can create file named EProjSimplex Dykstras algorithm process.
and call the EProjSimplex
algorithm in each iteration of the
function Ix Σtl EProfSimplex (v, k)
if nargin 2 k 1;
end;
ft-1;
length
v0 v-mean (v) k/n;
vmin min(v));
if vmin <
f = 1;
lambda_m 0;
while abs (f) 100-10
v1 v0 lambdaam;
posidx - v1>0;
noas - sum (posidx);
g -npos;
i sum (v1 (posidx)) k;
lambda m lamoda 2/gi
ft-ft+1;
if ft 100
x max (v1,
break;
end;
end;
max (v1, 0);
elsc
x = v0;
end;
The function ProjeoTTangent is simple projection function onto the tangent space.
function
Pro ToTangett (X)
c=s12c (x, 2);
cnxs (c, ;);
end
When running the code, please create four separate MATLAB files for RIMfactory, Dykstras, EProjSimplex new, and ProjToTannent,, and place them in the manopt folder following this structure:
-anopt;;
manifelds:
----muutii
RIMfactory;
Dykstras;
EProjSimplex_neww
Pro ProoiTTTangett;
Then you can call the functions in the general way as per manopt.
RIM_mani fold RIMfactory (n, c, row, upper, lower): ; problem.M RIM RIM_mani fold; problem. cost = 8 (X) ...;
problem. egrad (X) ...P Euclidean gradient rim, info rim,-l steepestdescent (problem)
56

---

Published as a conference paper at ICLR 2026
Furthermore, we provide reference code for the dual gradient and Sinkhorn algorithms, which allow the Retraction operation to be performed other ways. Overall, we still recommend using Dykstras algorithm under the Euclidean inner product for descent along the geodesics of the RIM manifold.
function dual _raddieg (2, 1, u, max iter)
[n, cl size (
1 ; u(:);
nu ones (n, 1
omega = ones (c, 1);
rho ones (c, 1);
step_size size .0;;
for iter - ::max_ites _iter
term - nu ones (1, ones (n, 1) omega' ones (n, 1)
rho';
__cureett max(term, 0);
grad nu = F current ones (c, 1) ones (n, 1) grad_omega F_current' ones (n, 1) grad_rhe -F_current's cnes 1)+1;
nu nu + step_size grad_nu; onega omega step_size size grad_omegaj rhe rhe step_size grad_rho;
omega max (omega, 0);
rhe max (rho, 0);
end
term
(1, cnes (h, 1) omega' (n, rho';
max (term, 0);
end
function on p (X, a, 1, N)
rng (1)
[n, c] size (X);
X;
u. vec = ones (n, 1);
vec ones(c, (c, 11;
V_vec ones (c, 1)
for 1 1:N ve = a // (K ve_ v_vec)
sum sum___t
(u_vec K), 1)'; vec max (1) ./ sum_P_t, ones 1));
sum sum_P_ P_t = sum vec K) vee' 1) v_vec min(u(:) ./ sum P_t, ones 1]];
F diag (u_
* diagiq_vec * v_vec);
B liehe = p oones (n, 1);
if norm (ones 1), 'fro') le-2 (1] a11/P_1cchosii-1o-22 all 21hhpllchh) break;
end
end
end
An alternative version of this paper is available on arXiv (Yuan et al., 2025).
57


Reviews:
[{"review_id": "cRy9YyORVl", "reviewer": "Reviewer_Y5CU", "summary": "A new relaxation is introduced for the indicator matrix optimization problem, leading to a manifold with a simple structure. Algorithms are designed based on this manifold.", "questions": "It is claimed after (1) that when $l=u=r$, the relaxation becomes $\\\\{X | X1_c=1_, X^T1_n=r, X>0\\\\}$. However, it seems to me that when $l=u=r$, the set in (1) is empty. What the tangent space of the manifold when $l=u$? Can you verify that the proposed retraction is still valid for this case?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "New retractions proposed and a class of efficient manifold algorithms are developed. Extensive numerical experiments have been conducted to validate the effectiveness of the proposed algorithms.", "weaknesses": "Presentation can be improved, e.g., \"Our Code is presented in Appendix H\", \"The proof is included in A.5\", etc.", "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": "tL9B9ZwA2R", "reviewer": "Reviewer_dz1d", "summary": "The paper introduces the Relaxed Indicator Matrix Manifold (RIM). The authors prove (M) is an embedded submanifold. They equip (M) with the Euclidean metric restricted to the manifold and derive a simple projection formula for the Riemannian gradient. \nFor retraction, they give two options: (1) a norm‑minimizing projection that they show is a geodesic for small steps (Dykstra‑style algorithm, Theorem 5, and (2) a Sinkhorn‑style mapping characterized by diagonal scalings.\n\nA central claim is lower per‑step complexity on RIM vs the doubly stochastic manifold (DSM) (Table 1).\n  \nExperiments: \n(i) retraction timing favors Dykstra as size grows (Table 2). \n(ii) “Experiment 2” compares RIM vs DSM on two problems: a convex norm‑approximation and TV image denoising.   \n(iii) “Experiment 3” applies RIM to Ratio‑Cut with closed‑form Euclidean gradient/Hessian (Appendix A.9) and shows losses/times vs several baselines\n(iv) “Experiment 4” reports clustering metrics (ACC/NMI/ARI) across datasets (Table 5). \n\nAuthors also provide a convergence proof of RIM‑RGD to stationarity at (O(1/T)) under standard smoothness and Armijo/Wolfe conditions.", "questions": "1. Can you provide any approximation or rounding guarantee that connects a stationary point on RIM to a discrete indicator solution (or to DSM) with a bounded increase in the objective?\n2. How should ([l,u]) be picked in practice without prior labels Do you recommend a data‑driven estimator beyond K‑means, and how sensitive are results across datasets?\n3. In Experiment 2, did both manifolds use identical stopping rules and the same retraction and line‑search settings If not, how do the results change when matched?\n4. Can authors provide concrete denoising metrics such as PSNR/SSIM/LPIPS/MSE for RIM vs DSM denoising experiments?", "limitations": "", "rating": 4, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "1. The RIM construction interpolates between single‑stochastic and doubly‑stochastic models via ([l,u]) bounds, letting practitioners encode prior class information. \n2. Clean Riemannian operators on RIM reduce to simple column‑mean corrections for the gradient and avoid DSM’s pseudoinverses, yielding the (O(nc)) vs (O(n^3)) complexity gap (Table 1).  \n3. Multiple retractions are implemented and compared; Dykstra is both fast at large sizes and, by Theorem 5, induces a geodesic in the small‑step regime.  \n4. Useful coverage of Ratio‑Cut with explicit Euclidean gradient and Hessian (Appendix A.9) that can be projected to RIM. \n5. Convergence guarantees for RIM‑RGD with Armijo/Wolfe steps are provided.", "weaknesses": "1. No approximation guarantees to the discrete indicator problem. The paper proves geometry and algorithmic convergence, not approximation quality: there is no integrality gap, rounding guarantee, or conditions under which RIM recovers the discrete optimum. The main theory sections and appendices focus on the toolbox, complexity, convergence, not on approximation bounds. A relaxation cannot just be fast, it must actually be a good approximation to the original problem (discrete indicator).\n2. TV denoising evidence is weak. The claim relies on an ~10% objective gap (1.05e5 vs 1.17e5) and “zoomed‑in” visual inspection; no PSNR/SSIM/LPIPS or even data‑term MSE are reported near those figures/tables [p.8, text around Table 3; Fig.4]. \n3. Compute fairness needs clarification. Speedups depend on stopping rules, retraction choice, and any hyperparameter tuning (e.g., (\\xi) in TV); these details are not fully standardized across RIM vs DSM in the description of Experiment 2. \n4. Choice of ([l,u]) is heuristic and dataset‑dependent. Appendix D.3.5 admits the difficulty, proposes using (n/c) or K‑means proportions, and shows a sensitivity table only for one dataset (MnistData05). \n5. Claims that “RIM images are clearer (regardless of noise level)” (page 47) are again qualitative and not well justified.", "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": "myrIbHzVUu", "reviewer": "Reviewer_m9ij", "summary": "The paper introduces the RIM (Relaxed Indicator Matrix) manifold as a generalization of the doubly stochastic and single stochastic manifolds. It proves Riemannian geometric properties of this set and develops a full Riemannian toolbox for optimization, including gradient, Hessian, and multiple retraction methods. The method is applied to large-scale tasks like Ratio Cut clustering and image denoising.", "questions": "1. How do you guarantee the solution remains in the interior of the RIM manifold? Do you clip or project back when the optimizer moves outside?\n\n2. Can the projection in Theorem 5 fail if the retracted point lies on the boundary (e.g., due to positivity or sum constraints)?", "limitations": "", "rating": 4, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "1. This paper proposes a manifold that interpolates between existing relaxations of indicator matrices.\n\n2. The authors develops practical Riemannian algorithms (including a retraction via projection).\n\n3. The authors shows speedups over doubly stochastic manifold optimization in experiments.", "weaknesses": "1. The paper defines the RIM manifold as a subset of Euclidean space with strict inequalities: M = {X | X1 = 1, l < Xᵀ1 < u, X > 0}\nHowever, this is an open set, and hence optimization problems posed over this set may have no solution, since the infimum may lie on the boundary (e.g., when X^T 1=u). The paper does not discuss this issue, and Theorem 1 claims it's an \"embedded submanifold\" without clarifying what happens near the boundary. This has major implications — e.g., the projection-based retraction (Theorem 5) may not be well-defined if the minimum lies at the boundary where feasibility breaks.\n\n2. The paper claims that optimization on the Stiefel manifold has time complexity \nO(n^3), but this is misleading. For typical tall matrices (n≫c), the complexity is only O(nc^2), as shown in Wen and Yin (2013). Please correct this and cite:\n\n       -Wen and Yin, A feasible method for optimization with orthogonality constraints, Math. Program. (2013)\n\nAlso, for retraction definitions, please cite standard sources:\n\n       -Boumal, “An Introduction to Optimization on Smooth Manifolds”, 2023, or\n\n       -Absil et al., “Optimization Algorithms on Matrix Manifolds”, Princeton, 2008\n\nIn Theorem 5, the projection used as a retraction may not always be well-defined — the feasibility set is open, and projection may fail to produce a point inside. The paper should acknowledge this and cite related projection-based retractions such as:\n\n       Absil & Malick, Projection-like retractions on matrix manifolds, SIAM J. Optimization, 2012.\n\n3. Every minor step is called a \"Theorem\", including obvious projections and inner product definitions. This makes the paper harder to read. Please distinguish between core theoretical contributions and auxiliary 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": {}}]

Author rebuttal:

Decision: Accept (Poster)