A Block Coordinate Descent Method for Nonsmooth Composite Optimization under Orthogonality Constraints
Authors: Ganzhao Yuan
Venue: ICLR 2026 Poster (2026)
Keywords: Orthogonality Constraints; Nonconvex Optimization; Nonsmooth Composite Optimization; Block Coordinate Descent; Convergence Analysis
Primary area: optimization

Abstract: Nonsmooth composite optimization with orthogonality constraints has a wide range of applications in statistical learning and data science. However, this problem is challenging due to its nonsmooth objective and computationally expensive, non-convex constraints. In this paper, we propose a new approach called \textbf{OBCD}, which leverages Block Coordinate Descent to address these challenges. \textbf{OBCD} is a feasible method with a small computational footprint. In each iteration, it updates $k$ rows of the solution matrix, where $k \geq 2$, by globally solving a small nonsmooth optimization problem under orthogonality constraints. We prove that the limiting points of \textbf{OBCD}, referred to as (global) block-$k$ stationary points, offer stronger optimality than standard critical points. Furthermore, we show that \textbf{OBCD} converges to $\epsilon$-block-$k$ stationary points with an iteration complexity of $\mathcal{O}(1/\epsilon)$. Additionally, under the Kurdyka-Lojasiewicz (KL) inequality, we establish the non-ergodic convergence rate of \textbf{OBCD}. We also demonstrate how novel breakpoint search methods can be used to solve the subproblem in \textbf{OBCD}. Empirical results show that our approach consistently outperforms existing methods.

Published as a conference paper at ICLR 2026
A BLOCK COORDINATE DESCENT METHOD FOR NONSMOOTH COMPOSITE OPTIMIZATION UNDER OR- THOGONALITY CONSTRAINTS
Ganzhao Yuan
Shenzhen University of Advanced Technology (SUAT), China yuanganzhao@ com
ABSTRACT
Nonsmooth composite optimization with orthogonality constraints has wide range applications in statistical learning and data science. However, this prob- lem is challenging due to its nonsmooth objective and computationally expensive, non-convex constraints. In this paper, we propose a new approach called OBCD, which leverages Block Coordinate Descent to address these challenges. OBCD is a feasible method with a small computational footprint. In each iteration, it updates & rows of the solution matrix, where k 2, by globally solving small nonsmooth optimization problem under orthogonality constraints. We prove that the limiting points of OBCD, referred to as (global) block-k stationary points, of- fer stronger optimality than standard critical points. Furthermore, we show that OBCD converges e-block-k stationary points with an iteration complexity of ((1/e). Additionally, under the Kurdyka-Lojj Lojassowawiwzz (KL) inequality, we estab- lish the non-ergodic convergence rate of OBCD. We also demonstrate how novel breakpoint search methods can be used to solve the subproblem in OBCD. Empir- ical results show that our approach consistently outperforms existing methods!"
INTRODUCTION
We consider the following nonsmooth composite optimization problem under orthogonality con- straints means define):
F(X) A f(X) + h(X), s.t. s.t. = XER
min
(1)
Here, n 2 r. n 2 2, and I, is a x r identity matrix. We do not assume convexity of f(X) and hX)). For brevity, the orthogonality constraints T I, in Problem (1) is rewritten as X St A (X Rarr
x I, where M A St r) is the Stiefel manifold in the literature (Edelman et al. 1998 Absil et al. 2008 Wen & Yin 2013 Hu et al. 2020). We impose the following assumptions
Problem (I) throughout this paper. (Asm-i) For any X and X where X and X only differ at most by k rows with & 2 2, we assume R2xr + IR is differentiable and H smooth with H R norr such that:
S(X+) S Q(
f(X) + (X+ X, VX(X)) + + X 1P
(2)
where || S Ly for some constant L, 0 and || A veo(X) Hvec X) Here, || so is the spectral norm of H. Notably, when H Ly Inr. this condition simplifies to the standard syysmoothssss (Nesterov) 2003). (Asm-ii) The function h(X) R + R is proper, lower semicontinuous, and potentially non-smooth. Additionally. it is coordinate-wist separable, such that l(() Ex, I((),)). Typical examples of ((X) include the l, norm l(X) = | with e
{0, 1), the capped-l; function h(X)= = Σ,, max([X.; with - > 0, and the indicator function
FFuture versions of this paper can be found at *Consider f(X)= = jur(XT -ur(XTCCD) CXD) ]| where D @ C, and C E RT*.
0364
E g are
symmetric. Clearly, f(X) satisfies (2) with equality, i.c., f(x+) Q X) for all X and x

---

Published as a conference paper at ICLR 2026
for non-negativity constraints h(X) (20(X). (Asm-iii) The following small-sized subproblem can be solved exactly and efficiently:
min VESI(R.R)
P(V) - V + +)
(3)
for any given Z E Rxxr P R and Q xx Here, we employ notational simplification by defining hVZ) A Ex, [1)) given the coordinate-wise separability of h(-). This assump- tion is analogous to the "prox- "rrondfdeendy"" condition in (variable-metric proximal gradient methods (Beck Teboulle, 2009) Raguet et al., 2013), but instead of evaluating standard proximal oper- ator for a single nonsmooth term in the full space, our subproblem jointly handles two nonsmooth components (the function h(-) and the orthogonality constraint) in a low-dimensional & x space.
Problem 1 is an optimization framework that plays crucial role in variety of statistical learning and data science models, such as sparse Principal Component Analysis (PCA) (Journée et al. 2010
Shalit & Chechik 2014), nonnegative PCA (Zass & Shashua 2006; Qian et al. [2021). deep neural
networks (Cogswell et al. 2016; Cho & LLee 2017; Massart Abrol, 2022] Huang & Gao 2023),
electronic structure calculation (Zhang al.
2014 Liu et al. 2014). Fourier transforms approx-
imation (Frerix & Bruna 2019), orthogonal nonnegative matrix factorization (Jiang et al. 2022
K-indicators clustering Diang al. 2016), and dictionary learning (Zhai et al. 2020).
1.1 MOTIVATING APPLICATIONS
Many machine learning and data science models can be cast as instances of Problem ① Below, we present two representative examples: In-regularized sparse PCA and Tn-reuularized sparse PCA. An additional example on nonnegative PCA is provided in Appendix Section G.]
Lo- Ro-Regularized Sparse PCA. Lo-regularized Sparse PCA (SPCA) is method that uses norm to produce modified principal components with sparse loadings, which helps reduce model com-
plexity and increase model interpretability (d d'Aspremont al., 2008 Chen al. 2016). It can be
formulated as: minx (t) - (X, CX) + X| where C = AA'A Ruxx is the covariance of the data matrix A E TR and À > 0.
L1 Regularrdd Sparse PCA. As the L1 norm provides the tightest convex relaxation for the Lo Lornorm over the unit ball in the sense of Lo- some researchers replace the non-convex and discontinuous Lo norm function with a convex but non-smooth function (Chen et al., 2016 Vu et al. 2013 Lu & Zhang, 2012). This leads to the following optimization problem of --grluulari SPCA: millees((n.))
CX) + X||X|11, where C Raxn is the covariance matrix of the data,
and A > 0.
1.2 RELATED WORK
We now present some related algorithms from the literature.
Minimizing Smooth Functions under Orthogonality Constraints. One of the main challenges in solving Problem ( stems from the nonconvexity of the orthogonality constraints. Existing ap- proaches for addressing this difficulty can be broadly grouped into four classes: (i) Geodesic-like methods (Abrudan et al.) 2008 Edelman et al., 1998 Absil et al. 2008 Computing exact geodesics typically involves solving ordinary differential equations, which can be computationally expensive. To avoid this, geodesic-like methods approximate the geodesic path by computing the geodesic log- arithm using simpler linear algebraic operations. (ii) Projection-like methods (Absil et al.) 2008 Golub & Van Loan 2013 Tiang & Dail (201)). These include techniques such as projection onto the nearest orthogonal matrix. polar decomposition. and QR-based projection. At each iteration, these methods descend along the Euclidean or Riemannian gradient direction and subsequently apply projection step enforce orthogonality. (iii) Multiplier correction methods (Gao et al.) 2018 2019, Xiao et al., 2022). These methods exploit the fact that the Lagrange multiplier associated with the orthogonality constraint is symmetric and admits a closed-form expression at first-order stationarity. They update the multiplier after achieving sufficient decrease in the objective, resulting in efficient
feasible or infeasible first-order methods. (iv) Landing methods (Ablin & Peyré 2022 Vary et al.) Ablin et al., 2024). These methods avoid explicit retractions by working in the ambient Euclidean space while adding a penalty that attracts iterates toward the orthogonal manifold. Each update com- bines descent direction for the objective with corrective term that reduces constraint violation,
2

---

Published as a conference paper at ICLR 2026
and, with appropriate step sizes, the iterates converge to points that are nearly orthogonal and nearly stationary for the original problem.
Minimizing Nonsmooth Functions under Orthogonality Constraints. Another major chal- lenge in solving Problem T arises from the nonsmoothness the objective function. Existing
approaches for handling this issue can be broadly categorized into four classes: (i) Subgradient methods (Hwang et al.] 2015 Li et al. [ et al.) 2024). These methods generalize gradi- ent descent to nonsmooth settings. Many of the previously mentioned geodesic-like and projection- based strategies can be incorporated into subgradient frameworks on manifolds. (ii) Proximal gra- dient methods (Chen et al., 2020 Li et al.) 2024b Lyu & Li 2025). These methods compute a descent direction by solving strongly convex subproblem over the tangent space, often using a semi-smooth Newton method. The resulting point then mapped back onto the manifold via a retraction to preserve orthogonality. (iii) Block Majorization Minimization (BMM) on Riemannian manifolds (Li et a 2024b; 2023] Breloy et al. 2021 Gutman Ho-Nguyen 2023). This class of methods iteratively constructs a tangential majorizing surrogate for a block of the objective, takes an approximate descent step in the corresponding tangent space, and retracts the iterate back to the manifold. (iv) Operator splitting methods (Lai & Osher 2014, Chen et al.] 2016: Zhang al.) 2019). These methods reformulate the original problem by introducing auxiliary variables and lin- constraints, decomposing it into simpler subproblems that can be solved separately and often exactly. Prominent examples include the Alternating Direction Method of Multipliers (ADMM) (He Yuan 2012), Riemannian ADMM (RADMM) (Li et al. 2024a), and Penalty- Splitting Method (PSM) (Yuan 2024 Chen 2012).
Block Coordinate Descent Methods. (Block) coordinate descent is a classical and powerful algorithm that solves optimization problems by iteratively performing minimization along (block) coordinate directions (Tseng Yun 2009: Xu & Yin 2013). The BCD methods have recently gained attention in solving nonconvex optimization problems, including sparse optimization (Yuan,
2024 Yuan et al. 2020), k-means clustering (Nie et al., 2022), structured nonconvex minimization Yuan 2023b a). recurrent neural network (Massart & Abrol 2022). and multi-layer convolutional
networks (Bibi et al., 2019 Zeng et al. 2019). BCD methods have also been used in (Shalit & Chechik 2014 Massart & Abrol 2022 for solving optimization problems with orthogonal group constraints. However, their column- holumn-wise BCD methods are Limited to solving smooth minimization problems with & = 2 and
n (refer to Section 4.2 in (Shalit & Chechik 2014)). Our row-wise BCD methods can solve coordinate- wise nonsmooth problems with /: 2 and r S n. The work of (Gao et al., 2019) proposes a parallelizable column-wise BCD scheme for solving the subproblems of their proximal linearized augmented Lagrangian algorithm. Impressive parallel scalability in a parallel environment of their algorithm is demonstrated. We stress that our row-wise BCD methods differ from the two column-wise counterparts.
Summary. Existing methods typically suffer from one or more of the following limitations: (i) they rely on full gradient information, incurring high computational costs per iteration; (ii) they do not accommodate coordinate-we nonsmooth composite objectives; (ii) they lack true descent properties and are often infeasible methods that only attain feasibility only in the limit; (ir) they often lack rigorous last-iterate convergence guarantees; (v) they provide only weak optimality results at critical points. * In contrast, our methods overcome these limitations by using tailored block coordinate descent framework for efficient composite optimization on the Stiefel manifold, with strong optimality and convergence guarantees.
1.3 CONTRIBUTIONS AND NOTATIONS
This paper makes the following contributions. (i) Algorithmically: We propose a Block Coordi- nate Descent (BCD) algorithm tailored for nonsmooth composite optimization under orthogonality constraints (Section 2. (ii) Theoretically: We provide comprehensive optimality and convergence analyses of our methods (Sections 3 and 4.. (ii) Empirically: Extensive experiments demonstrate that our methods surpass existing sol solutions in terms of accuracy and/or efficiency (Section 55).
We define [n] A (1.2.....n) m). and denote the Stiefel manifold as M 1 St(n, r). Matlab-style colon notation is used to describe submatrices. For matrix X E Raxr, let vec(X) Rarx1 denote the vector formed by stacking its columns, and let mat(x) E Rxr denote the inverse operator, such that mat x)) = X. We denote X 60 Y as the Kronecker product of the matrices X and Y. We
3

---

Published as a conference paper at ICLR 2026
use A + B and A IB to denote standard Minkowski addition and subtraction between sets A and B, and A @ IB and & è B to denote element-wise addition and subtraction, respectively. Additional notations are summarized in Appendix A.1
THE PROPOSED OBCD ALGORITHM
In this section, we introduce OBCD, Block Coordinate Descent algorithm for solving coordinate-
wise nonsmooth composite problems under Orthogonality constraints, as defined in Problem (
We start by presenting new update scheme designed to maintain the orthogonality constraint.
A New Constraint-Preservii Update Scheme. For any partition of the index vector [1, 21. n] into [B with B E No, B E N we define U3 E Ra and Up= E 1 as: (U2),z 1. = (Ugg)); jrr = { 1. = Therefore, we have the following variable s- 0, else.
0, else.
plitting for any X E RTXT; X I, (U,Uī + Un- ) (B,() + Un= ;)), where X(B.)) = U U E Rkxrr and x(B";))= Cl.X (n=-)jxr
In each iteration I, the indices {1, (1. of the rows of decision variable X E (n, r) are sepa- rated to two sets and Bº, where the working set with [B| and B = {1,2, ... n) B. To simplify notation,
B instead of B", as t can be inferred from the context. We only update k
rows of the variable X via X (B.:) = VX' ) for some appropriate matrix V E Ra*E The
following equivalent expressions hold:
X (B,:) (B;)
X (U& VU& U2- )
(4)
X = X U.( I))
(5)
We consider the following minimization procedure to iteratively solve Problem s.l. s. (V) E St(n,r) r), where X (V)
min V
X + UR(V Ix)UXX
(6)
The following lemma shows that the orthogonality constraint for X X+U - L.) can be preserved by choosing suitable V and X.
Lemma 2.1. (Proof in Appendix D. l We let B E [B. where the set (Bi, B2, ..., Boy } denotes
all possible combinations of the index vectors choosing k items from without repetition. We let
V E St(k, k). We define X+ & (V) & X + U.(V U (a) For any X E RRar, we have
[X [X+ x = x (b) X St(n,r), then X E St(n,r).
Thanks to Lemma 2.1 can now explore the following alternative formulation for Problem (6)
V E arg
arg min F(X (V)), s.t. V E St(k, k).
V
Then the solution matrix is updated via: X X (V*).
The following lemma offers important properties for the update rule X = X U.(V I Lemma 2.2. in
(.2) We define x X+U XVUU-V--IJUJ I. For any X E St(nn,))
V E St(k, k), B E J.h and symmetric matrix H E Rarrxr we have:
(a) |X X 2 I
L 122 where Q a (ZJ U²) HH Un), and Z A U.X Rxr. V. - X I
(b) X
X
2[I
(c) X X
IV
The Main Algorithm. The proposed algorithm OBCD is an iterative procedure that sequentially minimizes the objective function along block coordinate directions within a uub-manifold of M.
Starting with an initial feasible solution, OBCD iteratively determines working set B' using specif- ic strategies. then solves the small-sized subproblem in Problem (7) through successive Majoriza- tion Minimization (MM). This method iteratively constructs surrogate function that majorizes the objective function, driving it to decrease as expected (Mairal) 2013) Razaviyayn et al. 2013 Sun
et al. 2016 Breloy et al. 2021), and it has proven effective for minimizing complex functions.
4

---

Published as a conference paper at ICLR 2026
We now demonstrate how to derive the majorization function for (V)) in Problem 7 Initially,
for any X E St(n, r) and V E kt),,,),, we establish following inequalities: f(*" - f(x))
(V) X V
x X @ (U&(V - I v)(X))) + 2 V 1s ll1
XV 1. f X ) |ns)
+ / I 1 where step ① uses Inequality : step uses
Lemma 2.2 a); step @ uses 0 and Q N Q which can be ensured by choosing Q using one of
the following methods:
Q Q z Uz) H @ Us),
(8)
Q sI, with Q S S
where Z A X Then, we apply the MM technique to the smooth function f(X), while keep-
ing the nonsmooth component h(X) unchanged, leading to a function K(V; X', B) that majorizes F(2g = f(2 (V)):
(())(V)) S f(x) (V Is (XXX) bE) I Iella A(VU,X))
$ V Is -l (V. (X) (as)) I(VU
(10)
x(V:X'BB
where = f(X)) + X( - (X) ss) is a constant. Here, we use the
coordinate- separable property of h(-) follows: ((Ag(V) (Uge U. + VU!X)) h(U?. + xV) We minimize the upper bound of the right-hand side of Inequality (10) resulting in the minimization problem that V E arg minvest(k.)) K(V; X B), which can be effi- ciently and exactly solved due to our assumption.
Two simple strategies to find the working set B with |B| k can be considered. Random strategy: B is randomly selected from [BI, B2 .... Bes with equal probability 1/O,, (ii) Cyclic strategy: B' takes all possible combinations in cyclic order, such as B, - B2 - ... Bcg B₁
The proposed OBCD algorithm is summarized in Algorithm ! Importantly, OBCD is partial gra- dient method with low iterative computational complexity as only assesses rows of the Euclidean gradient of V VX(*) and the solution X to compute the linear term (IX (X')) V)
([V DI. X JRLsy V), as shown in Equation (10). Appendix C.3 details detais the complexity compari- son between OBCD and full gradient methods for some quadratic function f(X).
Algorithm OBCD: Block Coordinate Descent for Problem
1: Input: proximal parameter a > 0, initial feasible point X block size 2 2, =
2: for t = 0 to T do
3: (S1) Select a working set B E (1......n)*
n)* Denote B = 3² for simplicity.
4: (S2) Construct the matrix Q E " using 8 9
5: (S3) Define K(.,..-) as in Equation [T). and compute V as the global minimizer:
V arg miny K(V; B).
(11)
(Alternatively, find a local solution V such that IC (V. x B) < X(Lg;(() 6: (S4) ( V (B,:)
7:
end for
Solving the General OBCD Subproblems. The following lemma outlines key properties of the OBCD subproblems.
Lemma 2.3. (Proof in Appendix D.3) We define
= X and P A (X((X)(X*)]]
(at(QQee(Iee). - alg- We have:
(a) The subproblem in Equation 11 reduces to Problem with Q
(b) Assume that Formula ( is used to choose Q Problem (3) further reduces to the following problem: Vt E arg minvest((kk)) P(V) & (V,P) h(VZ). In particular, when h(X) & 0,
we obtain: V = (P))).. Here, PM (P) is the nearest orthogonality matrix to P.
5

---

Published as a conference paper at ICLR 2026
b). when k > 2, h(X)= O, and Q is chosen to be a diagonal matrix as in Equation (9), the subproblem V E arg min K B) in Algorithm | can be solved exactly and efficiently due to assumption, see Remark 2.6 (b) For general and 1(()), the subproblem may not admit a global solution. However, if a local stationary solution V satisfying K(V';; X B) X', B) can be found, then the sufficient descent condition remains valid, and convergence to a weaker optimality condition for the final solution X is still achievable (see Inequalities 44 (44)
Remark 2.4. (a) By Lemmna 2.3
Smallest Possible Subproblems When h = 2. We now discuss how to solve the subproblems exactly when ls The following lemma reveals an equivalent expression for any V St(2,
Lemma 2.5. (Proof in Appendix D.4) Any orthogonal matrix V E St(2, 2) can be expressed as cos(6) sin(()
V= V0th or V = Viel for some 0 E IR, where Viot - (- ( sin(O) cus(0) ) Vrel - cos(() sin(0)
siu(0) cas(0) )
We have del (Vrot) and det (Vre)) = for any 0.
Using Lemma
2.5 we can reformulate Problem 3 as the following one-dimensional problem:
Ó E arg min P(V), 8.L. E (V Vrof)
The optimal solution can be identified even h.)) * 0 using a novel breakpoint searching method, which is discussed later in Section B in the Appendix.
Remark 2.6. (i) Viott and Viel are called Givens rotation matrix and Jacobi reflection matrix respectively in the literature (Sun & Bischof 1995). Previous research only considered (Vist) for solving symmetric linear eigenvalue problems (Golub & Van Loan (013) and sparse PCA problems Shalit Chechik 2014), while we use {Vie V } for solving Problem (7). (ii) We show the necessity of using (o Voot
} in the following two examples of x optimization problems with orthogonality constraints: min F(V) - I|V AlE, and minvest(2.2 F(V) & +[V
Bl2 + +||[V where A = ( 21) and B = (12). The use of the reflection matrix Vg is essential in these examples because results in lower objective values. See Section C.1 in the Appendix for more details.
3 OPTIMALITY ANALYSIS
This section provides the optimality analysis for OBCD. First, we establish the completeness of the proposed update scheme, showing that OBCD can reach any feasible point from an arbitrary initialization. Second, we analyze the optimality conditions of both Problem ( and the associated subproblems of OBCD. Finally, by comparing these two sets of conditions, we derive hierarchy of optimality, illustrating how the algorithm's 's stationarity relates to that of Problem (
Basis Representation of Orthogonal Matrices. The following theorem shows that any orthog- onal matrix D E St(n,,)) n) and any point X E St(n,r) can be generated by composing simple d-imensioonl updates.
Theorem 3.1 (Basis Representation of Orthogonal Matrices), (Proof in Appendix 2. For all i E [C.], define W; & I, UB, (V, L)) Ug. V.U -UBUU
E.1 Assume k
where
V. St(k,)), Then:
(a) Any matrix D St(n,n) can be expressed as D We.,W, for suitable choice of W. (equivalently, V.). Furthermore, if'vi, V; = Ir, then D = In. (b) For any fixed reference point X E St(n, every X E rt(n,r) can be expressed as X We: ... W2W,X for suitable W..
The above representation for k = 2 can in fact be extended to any block size 2 2, as stated next.
Corollary 3.2. (Proof Appendix E.2) The conclusion of Theorem 3.7 extends to all k: 2
Remark 3.3. (i) We use both Givens rotation and Jacobi reflection matrices to compute D E St(n, n). This is necessary since a reflection matrix cannot be represented through sequence of rotations. (ii) The result of Corollary 3.2 indicates that the proposed update scheme X + X + UR(V -I, LT with V E St(k, k) as shown in Formula (5) can reach any orthogonal matrix X St(n,r) for any starting solution X St((,

---

Published as a conference paper at ICLR 2026
First-Order Optimality Conditions for Problem (1. We provide the first-order optimality
dition of Problem (T) (Wen & Yin 2013 Chen et al. 2020). We use (F) to denote the limiting
subdifferential of F(X (Mordukhovich 2006; Rockafellar & Wets. 2009), which is always non-
empty since F(X) is closed, proper, and lower semicontinuous. Given J(X) differentiable, we
have (x) = a( + h) (x) (h( (Rockafellar & Wets. 2009). We extend the def-
inition of limiting subdifferential to introduce FM (X) as the Riemanman limiting subdifferential
of F at X. defined as (MF(X) (F(X) @ (X(OF(X)) X), where @ is the element-
subtraction between sets.
Introducing a Lagrangian multiplier matrix A E Rx for the orthogonality constraint, we define the following Lagrangian function of Problem 1 f(X, A) = F(X) (I, X A). Notably, the matrix A is symmetric, as X is symmetric. We state the following definition of first-order optimality condition.
Definition 3.4. Critical Point Wen & Yin 2013 Chen et al. 2020). A solution X E St(n.r) r) is - (F( @ (X)FF(X) where (OF(X) O Moreover, the corresponding multiplier satisfies OF
critical point of Problem
( if: O 0M
X [ 1 X) (G X
A E [
Remark 3.5. The critical point condition in Lemma 3.4 can be equivalently expressed as Absil
et al. 2008: Jiang & Dai 2015 Liu et al. 2016): 0 E F (OF((X). Here, TxM is the tangent
space to M at X E M with TxM (Y gaxr X X 0}.
Optimality Conditions for the Subproblems. The Euclidean subdifferential of <C(V; X B²) (X)) where ((V
w.rl. V is given by GV) & Á(V) + ([(')
mat((Q + oIIg veo(v Ix and X = X U L X Using Lemma 3.4 we set the Riemannian subdifferential of K(V; W.t. V to zero and obtain the following first-order
optimality condition for V 0 E OMK((V;; is key ingredient in establishing the optimality hierarchy in Theorem 3.8 a) and the Riemannian subgradient lower bound in Lemma a).. a).
& G(V') e V'G(V') This inclusion
Optimality Conditions and Their Hierarchy. We introduce the following new optimality con- dition of block-k stationary points.
Definition 3.6. (Global) Block-k Stationary Point, abbreviated BSx-point. Let 0 and k 2 2. A solution X E St(n,r) is called a block-k stationary point if: VB E (8.) L arg minvest((k.)) K(V; B), where KC)::,)) is defined in Equation (10
Remark 3.7. BS.-poontt states that if we globally minimize the majorization function KC(; X, B), there is no possibility of improving the objective function value for K.V; X, across all (B.)
The following theorem establishes the relation between BSy-points, standard critical points, and global optimal points.
Theorem 3.8. (Proof in Appendix
E.3 We establish the following relationships:
(a) fcritical X} 2 (BS2-points X}.
(b) pBints X 2 Iglobaloptimal points X), where 32. (c) (BSie-points X} 2 [BS++1-points X), where k E (2.3...., n - 1}.
(d) The reverse of the above three inclusions may not always hold true.
Remark 3.9. (i) The optimality of BS2-points is stronger than that of standard critical points Wen Yin [013]] Chen et al. 2020; Absil et al. 2008). (ii Testing whether a solution X is a BSi-pooits deterministically requires solving all C₁ subproblems. However, by randomly selecting
the working set from the Cj possible combinations {B:) one can test whether X is -Sy-point in expectation.
7

---

Published as a conference paper at ICLR 2026
4 CONVERGENCE ANALYSIS
This section establishes the iteration complexity and non-crgodic (last-iterate) convergence rates of the proposed OBCD algorithm. We first prove sufficient descent property, followed by ergodic convergence rate typical in nonconvex optimization. We then analyze iteration complexity under the Riemannian subgradient condition, commonly used in nonsmooth manifold settings. Finally, derive a last-iterate convergence rate based on the KL inequality.
Throughout this section, assume that the working set is determined by random strategy and that the global minimizer V E arg minvesu((kk) (V;V; V; x','' 3²) can be computed. The algorithm OBCD then generates random output (V,X*+]) for t = 0, l...... 0o, depending on the realization of the random variable ç - (B?,B²,...,,B) We denote X as an arbitrary limit point of OBCD.
B²).
4.1 ITERATION COMPLEXITY
Initially, we introduce the notation of e-Sss-point as follows.
Definition 4.1. (c-BSt-point) Given any constant c point X is called an e-BSg-point if: C 2i dist(IR, arg miny XC X, B²))² S e, where K(+i.-) defined in Equation (10). Here, the set B2, ..., Boy }o denotes all possible combinations of the index vectors choosing k items from without repetition, and dist(E,E)) denotes the distance between two sets and E'.
Using the optimality measure from Definition 4.1 we establish the iteration complexity of OBCD. Theorem 4.2. (Proof in Appendix F.I We define ē a (F(x*) F(Xx)) 2 We have:
(a) The following sufficient decrease condition holds for all 1. 2 0:
+ 1 1 V' - L - F(X') - F(X*+).
(b) If the B² is selected from {B.).* randomly and uniformly, OBCD finds an e-BSa-point of Problem ④1 in at most T iterations in the sense of expectation, where T T [
Remark 4.3. Theorem 4.2 shows that OBCD converges to c-block-k: stationary points with iter-
ation complexity of ((1ee) which is typical for general nonconvex optimization.
Apart from Definition 4.1 another common optimality measure relies the Riemannian subgradi- ent. At the point V I. the Riemannian subdifferential of K(V;X', B²) is OMIC(I, X', BB) U. (De DB Us,, where D (XX(X++ @h( [ We next derive Riemannian subgradient lower bound in terms of the iterate gap.
Lemma 4.4. (Proof in Appendix F22 Riemannian Subgradient Lower Bound for the Iterates
Gap) Assume that F(.) is Cp-Lipschitz continuous on St(n, i.e., |G]| I CF for all X e St(n,r) r) and all E @F(X). We have:
(a) Egr-i [dist? (0, OMKK,, X (B))) 5o V I& where & - 72(C7. + +a²+I33)
(b) Ee (dist? (MF(()))) S Y Eg (dist?(0) 0M X B'))], where y 2 Ch. ////Ch=/ Ckk
Remark 4.5. The important class of nonsmooth & norm function h(X) || (Chen al. 2020 2024) satisfies the assumption made in Lemma 4.4
We establish the iteration complexity of OBCD using the optimality measure of Riemannian sub- gradient (Chen et al. 2020; Cheung et al. 2024 Li et al. (024b)).
Theorem 4.6. (Proof in Appendix F.3) We define č as in Theorem 4.2 and {ó, 7 in Lemma 4.4
OBCD finds an e-critical point of Problem I]. i.c., Bg (dist (0, F (xx+())))) I e, in at most
(0,
T ++1 iterations in expectation, where f [T] and T > 100
4.2 CONVERGENCE RATE UNDER KL INEQUALITY
We establish the non-ergodic convergence rate of OBCD using the a key tool in non-convex analysis (Attouch et al. 2010 Bolte et al. 2014
inequality, Liu et al. 2016).
Initially, we make the following additional assumption.
8

---

Published as a conference paper at ICLR 2026
Assumption 4.7. The function (. F(X) ((X)+MM)) (X) is Mardykka-ojasiewiizz (KL) function. Remark 4.8. Semi-algebraic functions constitute a broad class of KL. functions, including real polynomials, norm functions ||x||, with p 2 O, rank functions, and indicator functions of sets such as the Stiefel manifold and the positive semidefinite cone (Attouch et al. 2010))
We present the following useful proposition regrading to the KL function.
Proposition 4.9. (Kurdyka-Lojasieeiiee Property, see, e.g.. (....AAttocch al. 2010 Bolte et al. 2.14 Let F, mR - (-00, +o0] be a KI. function and X E dom Fi. Then there exist o E [), 1) E (0, +oo], neighhorhood Y of X and concave continuous function so()) ct- with e > 0 and E [0, ) such that for all X' E Y satisfying (X(X)) E (F,(X (Xx) + 1))
holds that dist(O, (x,()) vf (F((X')
(X 2 1.
Utilizing the HurdyyaaLLojaiieeee property, one can establish finite-length property of OBCD, result considerably stronger than that of Theorem [.2]
Theorem 4.10. (Proof in Appendix F:4 A Finite Length Property). We define E - (Eg [[[
I+]l)) 1/2 and D; Li=e Eijt1. Under the continuity assumption in Lemma 4.4 there exists sufficiently large t, such that, for all t,, we have
(a) It holds that (Ez-1)² < kpp: 9(+1), where f (F((() )) K 2/70 is
positive constant, y and Φ are defined in Lemma 4.4 and (p() is the desingularization function defined in Proposition 4.9,
(b) It holds that 2e= Es+1 S E, + 2kp... The sequence {E has the finite length property that D, A Ej= Fj+1 always upper-bounded by a certain constant for all / 2 l.
Finally, establish the last-iterate convergence rate for OBCD.
Theorem 4.11. (Proef in Appendix [.5). Based on the continuity assumption made in Lemma 4.4
for all t 2 Lu we have:
(a) If a = 0, then the sequence X' converges in a finite number of steps in expectation.
(b) o E (0. 2) then there exist > 0 and + E [0, 1) such that [g X S ctt.
(c) Ifo (2, 1), then there exist è > 0 such that Eg =|||''' [[]] & where 2 1-0 0.
2a-1
Remark 4.12. When F(X): is a semi-algebraic function and the desingularising function <p()) ct--a for some > 0 and a E 1), Theorem 4.17 shows that OBCD converges in finite iterations
when o with linear convergence when a E (0, 2). and sublinear convergence when o E (2,
for the gap IX' X= ||F in expectation. These results are consistent with those in (Attouch et al.
2010
5 EXPERIMENTS
This section presents numerical comparisons between OBCD and state-of-the-art methods on both real-world and synthetic data. We describe the application of To-regularized SPCA in the sequel, while additional applications for Li-reullarized SPCA and nonnegative PCA can be found in Ap- pendix Section G.2
Compared Methods on Legeegulad SPCA. We compare against three operator splitting methods: Linearized ADMM (LADMM) (Lai & Osher 2014, He & Yuan, 2012), Riemannian AD- MM (RADMM) (Li et al. 2024a). and the Penalty- based Splitting Method (PSM) YYuan 2024] 2012 Each method initialized with either random or identity matrix, yielding vari- LADMM(id), RADMM(dd), SPM(id), LADMM(rnd), RADMM (rnd), and PSM(()d). For OBCD, adopt random working strategy with identity initialization, denoted OBCD-
Chen
ants:
R(id).
Implementations. All methods are implemented in MATLAB on an Intel 2.6 GHz CPU with 32 GB RAM. However, our breakpoint searching procedure is developed in C++ and integrated into the MATLAB environment f as requires inefficient element-wing loops in native MATLAB. The
A lthough we prioritize accuracy over speed, the comparisons remain fair, as the other methods hased on matrix multiplication and SVD rely on highly optimized BLAS and LAPACK libraries.

---

Published as a conference paper at ICLR 2026
CADKINT
RADNIN
SPST
TRAINMM
RADION
CAOMIN
TE0T2.00 1000
156 0007
372 182
(30/997
559.776
2012
300 209-290
TD12 2012.000000 1000
1049 2003
0279 200
200
2140:271
12901225
6113 1000
10%
580
110 400
740
196651
159
6113 1000
2466
37 431
000
15
37.94 455
2-005-200000000 1003
299.113 113
1793999
250 120
1500777
E2305-2000 1003
1009
1241 523
1545-115
3393.85
2547
C Ferer 3099 1000
196 66.0
739 576
190 0.0
570 070
304
1639-055 657
196-97M
(2 CI 3993 1106
99% 06.0
000 400
10384 08.
63 6766
0009 467
1991981
48507T0
39:331
42199
2540.98
1000
19%
174 674
6715
198-971
1000
075
24 474
0700
676
19:9 470
vaada 900 1000
106
559 580
410 0000
5:0
199-077
rands 100% 1000
74.09
9600-077
1206
680
43:50426
RALLW0S
575
57R
LAURON
175
144
1451
2.40
17040
2040
100
10.43
1045
3700
1100 9113
7700 307
TIRR
10000
5200
TEDTD-90
5669 150
31
1000
$199.310
2599.218
237
TDTD
13391226
176
19989
553
A430 6403 1000
7799500
4929
2790716
30005500
43933513
1000060
A500 1000
13599 553
09991992
2550-199
MINIST 60000 754
1959 664
11:0 001
11296 715
1196 603
1804 653
0790 876
1996 1066
MNNST 00300 754
10994 655
55990 646
2619 804
1198 911
0000 035
1990
Clasts-
00-30000-100 1000
5469583
4399
11489/87
4399555
3892.583
0776233
1998-86
Clasti- 1000
13499 583
339991999
23066006
15495.080
19971978
19900000
Clee 1000
1958
1800
979
CTta: 1000 1000
15694
51964
1965-726
Table 1: Comparisons of objective values for Le-reuularized SPCA. The 1st, 2nd, and 3rd best results are colored with red, green and blue, respectively.
10
LADNDNUI RADRIME
LLANDAGI RADDOMOO)
LAUNDAII RAPAIMENO
Relitive Okjective
Debibe deyratee
LAUNDSOIS RATRIM
LAUNDAAIE RADDOM 1Nt
LAUNDAMOO
OACT-804
Time (seconds)
60
31
60
20
20 50 (secords)
(a) w. wlla22477-300
(b) TDT2 500-1000
(c) 20 8000--0000000000 1030
(d) sector 6412-100
Figure 1: The convergence curve for solving Lo-regularized SPCA with X = 100. No matter how long the algorithms run, the other methods remain trapped in poor local minima.
code for all three applications used reproduce the experiments can be found in the supplemental material.
Experiment Settings. We compare objective values F(X) for different methods after running for 30 seconds. For numerical stability in reporting the objectives, we use the count of elements with absolute values greater than threshold of 10-6 instead of the original E, norm function | We set a 10-5 for OBCD. Full-gradient methods have higher per-iteration complexity but require fewer iterations, while OBCD, partial-gradient method, has lower per-iteration costs but needs more iterations. Thus, we compare based on CPU time rather than iteration count.
Experiment Results. Table I and Figure I display accuracy and computational efficiency results for Lo-regularized PCA, yielding the following observations: (i) OBCD-R delivers the best per- formance. (ii) Unlike other methods where objectives fluctuate during iterations, OBCD-R mono- tonically decreases the objective function while maintaining the orthogonality constraint. This is because OBCD is a greedy descent method for this problem class. (iii) While other methods of- ten get stuck in poor local minima, OBCD-R escapes from such minima and generally finds lower objectives, aligning with our theory that our methods locate stronger stationary points.
6 CONCLUSIONS
In this paper, we introduced OBCD, a new block coordinate descent method for nonsmooth compos- ite optimization under orthogonality constraints. OBCD operates k rows of the solution matrix, offering lower computational complexity per iteration for 2 2. We also provide a novel optimality analysis, showing how OBCD exploits problem structure to escape bad local minima and find bet- ter stationary points than methods focused on critical points. Under the Kurdyka-Lojasiewin inequality, we establish strong limit-point convergence. Additionally, we show how novel break- point search methods can be used to solve the subproblem when k = 2. Extensive experiments demonstrate that OBCD outperforms existing methods.
(KL)
10

---

Published as a conference paper at ICLR 2026
REFERENCES
Pierre Ablin and Gabriel Peyré. Fast and accurate optimization on the orthogonal manifold without retraction. In International Conference on Artificial Intelligence and Statistics, pp. 5636-5657.
PMLR, 2022.
Pierre Ablin, Simon Vary, Bin Gao, and Pierre-Antoine Absil. Infeasible deterministic, stochastic, and variance-reduction algorithms for optimization under orthogonality constraints. Journal of Machine Learning Research, 25(389):1-38, 2024.
Traian E Abrudan, Jan Eriksson, and Visa Koivunen. Steepest descent algorithms for optimiza- tion under unitary matrix constraint. IEEE Transactions on Signal Processing, 56(3):1134-1)) 1147, 2008.
Pierre-Antoine Absil, Robert E. Mahony, and Rodolphe Sepulchre. Optimization Algorithms on Matrix Manifolds. Princeton University Press, 2008.
Hedy Attouch, Jérôme Bolte, Patrick Redont, and Antoine Soubeyran. Proximal alternating mini- mization and projection methods for nonconvex problems: An approach based on the kurdyka- lojasiewicz inequality. Mathematics of Operations Research, 35(2):438 457, 2010.
Amir Beck and Mare Teboulle. fast iterative dhrikkag algorithm for linear inverse problems. SIAM journal on imaging sciences, 2(1):183-02,, 2009.
Adel Bibi, Bermard Ghanem, Vladlen Koltun, and René Ranftl. Deep layers as stochastic solvers. In International Conference on Learning Representations (ICLR), 2019.
Jérôme Bolte, Shoham Sabach, and Marc Teboulle. Proximal alternating linearized minimization for nonconvex and nonsmooth problems. Mathematical Programming. (46(1-2):459 494, 2014.
Arnaud Breloy, Sandeep Kumar, Ying Sun, and Daniel P Palomar. Majorization-miniizaatiii minimization on the stiefel manifold with application to robust sparse pea. IEEE Transactions on Signal Processing. 69:1507-1520, 2021.
Shixiang Chen, Shiqian Ma, Anthony Man-Cho So, and Tong Zhang. Proximal gradient method for nonsmooth optimization over the stiefel manifold. SIAM Journal on Optimization, 30(1): 210-20 (00-2239.220 2020.
Shixiang Chen, Shiqian Ma, Anthony Man-Cho So. and Tong Zhang, Nonsmooth optimization over the stiefel manifold and beyond: Proximal gradient method and recent variants. SIAM Review, 66 (2):319-352, 2024.
Weigiang Chen, Hui Ji, and Yanfei You. An augmented lagrangian method for & e-rreullrrzed opti- mization problems with orthogonality constraints. SIAM Journal on Scientific Computing, 38(4); 2570-B592. 2016.
Xiaojun Chen. Smoothing methods for nonsmooth, nonconvex minimization. Mathematical pro- gramming, 134:71-99, 2012.
Andy Yat-Ming Cheung, Jinxin Wang. Man-Chung Yue, and Anthony Man-Cho So. Randomized submanifold subgradient method for optimization over stiefel manifolds. arXiv preprint arX- iv:2409.01770, 2024.
Minhyung Cho and Jaehyung Lee. Riemannian approach to batch normalization. Advances in Neural Information Processing Systems, 30, 2017.
Michael Cogswell, Faruk Ahmed, Ross B. Girshick, Larry Zitnick, and Dhruv Batra. Reducing overfitting in deep networks by decorrelating representations. In Yoshua Bengio and Yann LeCun (eds.), International Conference Learning Representations (ICLR), 2016.
Alexandre d' Aspremont, Francis Bach, and Laurent El Ghaoui. Optimal solutions for sparse princi- pal component analysis. Journal of Machine Learning Research, 9(7), 2008
Alan Edelman, Tomás A. Arias, and Steven Thomas Smith. The geometry of algorithms with orthog- onality constraints. SIAM Journal Matrix Analysis and Applications, 20(2):303-333, 1998.

---

Published as a conference paper at ICLR 2026
Thomas Frerix and Joan Bruna. Approximating orthogonal matrices with effective givens factoriza- tion. In Interaational Conference on Machine Learning (ICML), pp. 1993-2001, 2019.
Bin Gao, Xin Liu, Xiaojun Chen, and Ya-xiang Yuan. A new first-order algorithmic framework for optimization problems with orthogonality constraints. SIAM Journal on Optimization, 28(1): 302-332, 2018.
Bin Gao, Xin Liu, and Ya-xiang Yuan. Parallelizable algorithms for optimization problems with orthogonality constraints, SIAM Journal on Scientific Computing, (13)::1949-A198) 2019.
Gene H Golub and Charles F Van Loan. Matrix computations. 2013.
David H Gutman and Nam Ho-Nguyen. Coordinate descent without coordinates: Tangent subspace descent riemannian manifolds. Mathematics of Operations Research, 48(1):127-159, 2023.
Bingsheng He and Xiaoming Yuan. On the (1)/n) convergence rate of the douglas- alter-
nating direction method. SIAM Journal on Numerical Analysis, 5002):700-709, 2012.
Jiang Hu, Xin Liu, Zai-Wen Wen, and Ya-Xiang Yuan. brief introduction manifold optimiza- tion. Journal of the Operations Research Society of China, 5(2):199-244 2020.
Feihu Huang and Shangqian Gao. Gradient descent ascent for minimax problems riemannian manifolds. IEEE Trans. Pattern Anal. Mach. Intell., 45(7):8466-8476, 2023.
Scong Jae Hwang, Maxwell D. Collins, Sathya N. Ravi, Vamsi K. Ithapu, Nagesh Adluru, Sterling C. Johnson, and Vikas Singh. A projection free method for generalized eigenvalue problem with a nonsmooth regularizer. In IEEE International Conference on Computer Vision (ICCV), pp. 1841- 2015.
Bo Jiang and Yu-Hong Dai. A framework of constraint preserving update schemes for optimization on stiefel manifold. Mathematical Programming, 153(2):535-575, 2015.
Bo Jiang, Ya-Feng I. Liu, and Zaiwen Wen. f--norm regularization algorithms for optimization over permutation matrices. SIAM Journal on Optimization, 26(4):2284 2313, 2016.
Bo Jiang, Xiang Meng, Zaiwen Wen, and Xiaojun Chen. An exact penalty approach for optimization with nonnegative orthogonality constraints. Mathematical Programming, pp. 1-43, 2022.
Michel Journée, Yurii Nesterov, Peter Richtárik, and Rodolphe Sepulchre. Generalized power method for sparse principal component analysis. Journal of Machine Learning Research, 11
(2), 2010.
Rongjie Lai and Stanley Osher. A splitting method for orthogonality constrained problems. Journal of Scientific Computing, 58(2):331-445 2014.
Jiaxiang Li, Shiqian Ma, and Tejes Srivastava. A riemannian alternating direction method of multi- pliers. Mathematics of Operutions Research, 2024a.
Xiao Li, Shixiang Chen, Zengde Deng, Qing Qu. Zhihui Zhu, and Anthony Man-Cho So. Weakly convex optimization over stiefel manifold using riemannian subgradient-tyy methods. SIAM Journal on Optimization, 31(3):1605-16 2021.
Yuchen Li. Laura Balzano, Deanna Needell, and Hanbaek Lyu. Convergence and complexity of block mmaiorizttio-miinnnzzaoooo minimization for constrained arXiv:2312.103300 2023.
optimization. arXiv preprint
Yuchen Li, Laura Balzano, Deanna Needell, and Hanbaek Lyu. Convergence and complexity guar- antee for inexact first-order riemannian optimization algorithms. In International Conference on Machine Learning (ICML), 2024b.
Huikang Liu, Weijie Wu, and Anthony Man-Cho So. Quadratic optimization with orthogonality constraints: Explicit lojasiewicz exponent and linear convergence line-search methods. In International Conference Machine Learning (ICML), pp. 1158-1167, 2016.
12

---

Published as a conference paper at ICLR 2026
Xin Liu, Xiao Wang, Zaiwen Wen, and Yaxiang Yuan. On the convergence of the self-consistent density functional theory. SIAM Journal on Matrix Analysis and Applications, 35(2):546-33) 558, 2014.
field iteration in
Zhihua Lu and Yi Zhang. An augmented lagrangian approach for sparse principal component anal- 2012.
ysis. Mathematical Programming,
Hanbaek Lyu and Yuchen Li. Block intiorizatio-mmmmmmzzan with diminishing radius for con- strained nonsmooth nonconvex optimization. SIAM Journal on Optimization, 35(2):842-871, 2025.
Julien Mairal. Optimization with first-order surrogate functions. In International Conference on Machine Learning (ICML), volume 28, pp. 783-791, 2013.
Estelle Massart and Vinayak Abrol. Coordinate descent on the orthogonal group for recurrent neu- ral network training. In Proceedings of the AAAI Conference on Artificial Intelligence (AAAT), volume 36, pp. 7744-7751, 2022.
Boris S. Mordukhovich. Variational analysis and generalized differentiation i: Basic theory. Berlin Springer, 330, 2006.
Yurii Nesterov.
Introductory lectures on convex optimization: A basic course, volume 87. Springer Science & Business Media, 2003.
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-2388 2022.
Yitian Qian, Shaohua Pan, and Lianghai Xiao. Exact penalty methods for minimizing a smooth function over the nonnegative orthogonal set. arXiv, 11 2021.
Hugo Raguet, Jalal Fadili, and Gabriel Peyré. A generalized forward-backward splitting. SIAM Journal on Imaging Sciences, 6(3): (139:1199-11 199--1226,2,,,,, 2013.
Meisam Razaviyayn, Mingyi Hong, and Zhi-Quan Luo. A unified convergence analysis of block successive minimization methods for nonsmooth optimization. SIAM Journal on Optimization, 23(2):1126-1153 2013. 2013.
R. Tyrrell Rockafellar and Roger J-B. Wets. Variational analysis. Springer Science & Business Media, 317, 2009.
Uri Shalit and Gal Chechik. Coordinate-descent for learning orthogonal matrices through givens rotations. In International Conference on Machine Learning (ICML), pp. 548-556. PMLR, 2014.
Xiaobai Sun and Christian Bischof. A basis-kernel representation of orthogonal matrices. SIAM Journal on Matrix Analysis and Applications, 16(4):1184-1199 184-119,, 1995.
Ying Sun, Prabhu Babu, and Daniel P Palomar. Majorization-mininiin
algorithms in signal processing, communications, and machine learning. IEEE Transactions on Signal Processing, 65 (3):994-816, 2016.
Paul Tseng and Sangwoon Yun. A coordinate gradient descent method for nonsmooth separable minimization. Mathematical Programming, 117(L-2):877-433, 2009.
Simon Vary, Pierre Ablin, Bin Gao, and P.-A. Absil. Optimization without retraction on the random generalized stiefel manifold. In International Conference on Machine Learning (ICML).
Vincent Q Vu, Juhee Cho, Jing Lei, and Karl Rohe. Fantope projection and selection: A near-optimal convex relaxation of sparse pca. Advances in Neural Information Processing Systems (NeurIPS), 26, 2013.
Zaiwen Wen and Wotao Yin. A feasible method for optimization with orthogonality constraints. Mathematical PraiicaallPrrreree Programming, 142(1-2):397, 2013.
13

---

Published as a conference paper at ICLR 2026
WikiContributors.
Quartic equation: httpa://en. http:///aiiiHHiooaaa.oooosaa///
equation
pp. Last edited in March, 2023.
Nachuan Xiao, Xin Liu, and Ya-xiang Yuan. A class of smooth exact penalty function methods for optimization problems with orthogonality constraints. Optimization Methods and Software, 37 (4):1205-1241, 2022.
Yangyang Xu and Wotao Yin. A block coordinate descent method for regularized multiconvex optimization with applications nonnegative tensor factorization and completion. SIAM Journal on Imaging Sciences, 6(3):1758-1793) 2013.
Ganzhao Yuan. Coordinate descent methods for de minimization: Optimality conditions and global convergence. Proceedings of the AAAI Conference on Artificial Intelligence (AAAI), 2023a.
Ganzhao Yuan. Coordinate descent methods for fractional minimization. In International Confer- ence on Machine Learning (ICML), 2023b.
Ganzhao Yuan. Smoothing proximal gradient methods nonsmooth sparsity constrained optimiza- tion: Optimality conditions and global convergence. In International Conference Machine Learning, 2024.
Ganzhao Yuan, Li Shen, and Wei-Shi Zheng. A block decomposition algorithm for sparse opti- mization. In ACM SIGKDD Conference Knowledge Discovery and Data Mining (SIGKDD), 2020.
Ron Zass and Amnon Shashua. Nonnegative sparse PCA. Advances in Neural Information Process- ing Systems (NeurIPS), pp. 1561-1568, 2006.
Jinshan Zeng, Tim Tiz-Kit Lau, Shaobo Lin, and Yuan Yao. Global convergence of block coordinate descent in deep learning, In International Conference Machine Learning (ICML), pp. 7313- 7323. PMLR, 2019.
Yuexiang Zhai, Zitong Yang, Zhenyu Liao, John Wright, and Yi Ma. Complete dictionary learning via 14-norm maximization over the orthogonal group. Journal of Machine Learning Research, 21 (165):1-68, 2020.
Junyu Zhang, Shiqian Ma, and Shuzhong Zhang. Primal-dual optimization algorithms over rieman- nian manifolds: an iteration complexity analysis. Mathematical Programming, pp. 1-46, 2019.
Xin Zhang, Jinwei Zhu, Zaiwen Wen, and Aihui Zhou Gradient type optimization methods for electronic structure caleulations. SIAM Journal on Scientific Computing, 36(3):265-289, 2014.
LLM USAGE
large language model (LLM) was used to assist refining the writing of this paper.
ACKNOWLEDGMENTS
This work was supported by Guangdong Natural Science Funds for Distinguished Young Scholar (2018 88030306002).
14

---

Published as a conference paper at ICLR 2026
Appendix
The appendix section is organized follows.
Section
A
covers notations, technical preliminaries, and relevant lemmas.
Section
B
shows how solve the subproblem when k
Section
offers further discussions on the proposed algorithm.
Section
D
contains proofs from Section 2
Section contains proofs from Section
Section contains proofs from Section
Section presents additional experiment details and results.
A NOTATIONS, TECHNICAL PRELIMINARIES, AND RELEVANT LEMMAS
A.1 NOTATIONS
Throughout this paper, M A St(n, r) denotes the Stiefel manifold, which an embedded subman- ifold of the Euclidean space gaxr, Boldfaced lowercase letters denote vectors and uppercase letters denote real-valued matrices. We adopt the Matlab colon notation to denote indices that describe submatrices. For given natural numbers n and k use [B1, B2,.... Beg } to denote all the possi- ble combinations of the index vectors choosing items from n without repetition, where C the total number of such combinations and B, N Vi E [C,]. For any one-dimensional function p(t) R + I. define: p(#r v) A minfp(x y). p( y)}. We use the following notations in this paper.
[n]: 2,.,,..,ni
|x|: Euclidean norm: ||x|| = ||x|22 = xx
X₁ the i-th element of vector
X,, or Xij the (ith, jih element of matrix X
voc(X): (ec(X) Rerx1, the vector formed by stacking the column vectors of X
mat(x): mat (x) E Rexr, Convert x E Rarxl into matrix with mat(voo(X))
X the transpose of the matrix X
sigu(t): the signum function, sign(() = if and sign(() =
otherwise
det(D): Determinant of square matrix D E -
Ck the number of possible combinations choosing items from n without repetition matrix of size
Oz,r
the subscript is omitted sometimes
1 I, E IR" Identity matrix
X ž O(or > 0) the Matrix X is symmetrie positive semidefinite (or definite)
tr( A)) Sum of the elements on the main diagonal X: tr(A) Σ; A4,
X, Y)
Euclidean inner product, i.e., (X, Y) E. X, Y
X x Y Kronecker product of X and Y
| F Frobenius norm: (Ee X2
X app Operator/Spectral norm: the largest singular value of X
XX the number of non-zero elements in the matrix X
|X|]: the absolute sum of the elements in the matrix X with ||X||] Ei, |X,,yl
m the capped-é, norm of X with || max([X), 7) = Ei, max( (Xisl:7) 7)
V ( Euclidean gradient of f(X) at X
MM (X) Riemannian gradient of f(X) X
15

---

Published as a conference paper at ICLR 2026
(
limiting Euclidean subdifferential of F(X) at X
QMF(X)
liiiiiig Riemannian subdifferential of F(X) at X
3 x the indicator function of a set E with 1= if E E and otherwise +00
220 X
indicator function of non-negativity constraint with (20 (X) =
{
u.
2 0.
oc.
cbe.
Ps
Orthogonal projection of with P=(Z) arg minx= I X
PM Y Nearest orthogonal matrix of Y with I", (Y) arg migx |X Y
dist( [, =') the distance between two sets with dist EE _ inxxe=xx'EE | X
IB, A IB: standard Minkowski addition and subtraction between A and IB
A + B, A IB: element-win addition and subtraction between sets A and B
[[F((X)]E ||:: the distance from the origin to (F)) with [|OF(X)]IF inf ||]|F
A.2
TECHNICAL PRELIMINARIES
As the function F(.) can be non-convex and non-smooth, we introduce some tools in non-smooth analysis (Mordukhovich) 2006 Rockafellar & Wets., 2009). The domain of any extended real- valued function F R= + (-o0, +oo| is defined as dom (F) & [X gaxr [F(X) < +00).
The Fréchet subdifferential of F at X E dom(F)) is defined as
OF(X) A (E E Ruxr:
lim inf F(Z)-F(XX E2 X) 0). Z-X 2
-
while the limiting subdifferential of F(X) at X E dom(F) is denoted as
OF(X) - (EER" EX + X, F (X') - F(X), E OF((X) -, V...
We denote VF(X) the gradient of F(.) at X in the Euclidean space. We have the following relation between (F(X), (F(X), and (F( (i) holds that (F(x) OF(X). (ii) If the function F(.) is convex, OF(X) and (F(X) essentially the classical subdifferential for convex functions, i.e.,
OF(x) OF(X) (EER** F(()) 2 F(X) + (E. X), VZ E "
(iii) If the function F(-) is differentiable, then OF(X) aF = (VF()))
We need some prerequisite knowledge in optimization with orthogonality constraints (Absil et al.) 2008 The nearest orthogonality matrix to an arbitrary matrix E Rn is given by PM UV where Y = (Diag(s) VT is the singular value decomposition of Y. We use A to denote the limiting normal cone M X, leading
(x) (X) {Z E nax, (Z, YY
The tangent and norm space to M at X E M are denoted as Tx and NxM, respectively. For a given X E M, we let Ax (Y) T + Y for Y E Raxr and we have Tx. = = e Raxr (x
0} and Nx.M (2XA A A Rrxr). For any non-convex and non- smooth function F(X). we use FMF to denote the limiting Riemannian gradient of F(X) at
X. and obtain 0MF(() = (OFM(OF(X)), We denote OF(X) e X(OF(X)) A E
G XG G X, OF()).
A.3 RELEVANT LEMMAS
We offer a set of useful lemmas, each of which stands independently of context and specific method- ology.
Lemma A.1. Let k 2 2 and W E R"x", Of W. for all B E (B.) then W 0. Here, the set {B,, BB, ..., Bes represents all possible combinations of the index vectors choosing items from without repetition.
16

---

Published as a conference paper at ICLR 2026
Proof. The proof is straightforward and relies on elementary reasoning.
Notably, the conclusion of this lemma does not necessarily hold if [B] = k = This is because any matrix W E _x with W. = 0 for all i [ni] satisfies the condition of this lemma but is not necessary zero matrix.
Lemma A.2. For any matrices A R and E Rkxk we have:
||
F S 2 ClI + C CT
Proof. We derive: |A AT = (| (C CT) (C- ACT - ||C
C + | A F = 2 - ClI + |C - C where step D uses the triangle inequality.
A
Lemma A.3. Let E R, and A E g2x2 be any skey-symmettric matrix with A
-A. We have:
det ((I2 _A [I, [A]) =.
Proof. Since A is two-dimensional matrix, it can be expressed in the form: A a E R. Letting b = a. we derive:
0 0) for some
Q
L _ A) (I _ A) @ (1 (t 2
_ - (1
1-b² 26
1+62
25 1-
where step I uses _ A b): step @ uses the fact that (2
ad-be be (-c -b] -1 for all
b. c, d E R. We further obtain: det (Q) @ 1-1² 1-6²
2b -bb
(1-622445) (1+62)2
4b2
(1462y3 = 1,
(1+82)
1+52 1462 14t 1+
where step uses the fact that det( a )) ad be for all a, b.b. R.
Lemma A.4. For any W pa we have
ct Σ
+
EC Σ W
3.J
Here, the sei (Bi, B2,... ... Bog represents all possible combinations of the index vectors choosing items from without repetition.
Proof. For any matrix W E Rax", we define: A diag( W) E R and W W - -Diagww).
We have: W
Diag( w)
+ W', this leads to the following decomposition:
WU&, = Σ |U B, (Diag(w) w) + W 1F
B
UUB B. Diag( w) U + D
B,
W' 1..
(12)
T2
We first focus on the term r We have:
r = -= Ch IU Diag(w) w Us, 12 2 E |a.
C
|| C+
Σ. W (13)
where step ① uses the fact that ||B Diag( B [[Diag(w) w) ]es 11 | l22 for any B E (3
step @ uses the observation that w; appears in the term Σ - || 112 total of (C,, 3 times, which can be deduced using basic induction.
We now focus on the term l2. Noticing that W # = 0 for all i E [n]. we have:
r2 1= == JUT WU Σ.
E C 3
2 c
(14)
where step T uses the fact that the term E each including {i.}} with + j: step uses
W ,,))
| comprises Cx distinet patterns,
Σ Eiyz,
In view of Equalities
13 and 14 complete the proof of this lemma.
17

---

Published as a conference paper at ICLR 2026
Lemma A.5. Assume QR X R=x", where Q St(n, n) and R is a lower triangular matrix with R. = 0 all <. If E St(n,n). then R a diagonal matrix with R E {-1,+11)foo all i [n].
Proof. We derive: RR (QX) QxX I. where step 1 uses R = Q' X; step @ uses X E St(n, n) and Q St(n, n). First, given [|[ and R n) 0, we have R1,1 E +1]. Second. we have [| and R(1,:) ((:;,) 0, leading to R1.2 0 and R2,2 E [-1,+1}. Finally, using similar recursive strategy, we conclude that is diagonal matrix with R,,; E (-1,+1] +1] for all E [m).
Lemma A.6. We define TxM A (Y E R [Ax(Y) 0} and Ax & X + + For any E Raxx and X St(n,, k), we have:
(G 2 (XAx(G))= = arg min YETxM
|Y - G
Proof. The conclusion of this lemma can be found in Absil et al. 2008 For completeness, we present short proof.
Consider the convex problem: Y = arg miny IY - Gll s.t. X + Y 0. Introduc-
ing multiplier E R.xr for the linear constraints leads to the following Lagrangian function:
L(Y; A) = Y - G + (X + YT, We derive the subsequent first-order optimali-
ly condition: 2(Y G) + X + AT) O, and X + Y
Given A is symmet-
ric, we have Y = G XA. Incorporating this result into X Y O. we obtain:
( - XA) + (G-XA)
O. Given X St(n,r), r), we have X A + - A - leading to: A 1[XX G Therefore, the optimal solution Y can be computed as
O,
Y XA = --G-X(XXGGGG
Lemma A.7. Consider the following problem: minx F.(X) & F(X) + (M(X), where F(X) is defined in Equation 7 For any X St(n, r), holds that
dist(0, (.) I (M)())
Proof. We let e (F(X) and define (x A G
Recall that the following first-order optimality conditions are equivalent for all X E St(n,r): (0 (X)) e (0 E (FF)((F(X))). Therefore, we derive:
(0,
infyeoF((x) || infyer
|Y|F
e
())
|| [XA(
|G +
[| [XXT) -
|G - X
where step ① uses G E (X)::; step uses Lemma step uses the definition of Ax (G));; step @ uses the identity that G [X + G (G XG step @ uses the norm inequality and fact that the matrix I x only contains eigenvalues that are 1 or 1.
Lemma A.8. Assume cos(0) # O. Any pair of trigonometric functions (ins()),sin(0)) can be represented as follows:
a) cos(0) (o+ and sin(O)
tan(e)
/ VI+tan"
18

---

Published as a conference paper at ICLR 2026
b) cos(0) -1 and sin(0) (I+
tan(
(0)
Proof. For all values of & where cos(0) + 0. the trigonometric functions {sin cos( t0) are well-defined. Utilizing the identity sin? (0) + cos (0) and tan(0) cos(O) sin(0), we
derive: (tan(f) cos cos)(0) = 1. Consequently, we find: cos(O) v (0)+1
Finally, we
can express sin as sf) t@n cos(0) (8)+1
tan(0)
Lemma A.9. Assume (F4-1)2 S Falpe P1+1) and Pt 2 p++1, where {ee Pt ptli are two negative sequences. For all i 2 1, we have: Etsi Ex+1 S E. + 2p.
Proof. We define wi & Pt- P-+1. Weleiii<I.
First, for any 2 1, we have:
EL We
Et T Pc+1) ps -PT+1 Ph
(15)
where step uses p. 2 for all i. Second, we obtain:
M Ino.
E141
v V
(F.) (w)) (2nx,, Va: >
V Es + wi 1/(2a), Va > 0.
(16)
Here, step ① uses (E++1)² < Ex(p: Pc+1) and Pt Pi+1: step @ uses the fact that ab < şa² + b² 6² for all a > 0; step @ uses the fact that Va + b S Va vb for all a, 6 2 0.
Assume / > O. Telescoping Inequality 16 over t from i to T, have:
E. wt I1((a)
1ET Ec+1} / 1 Ed}
{ETT + Σ E++1} v / Ex + Σt=
11'-1
E
ETtI -V
- Es + / ET-1 E
- --VE + (1 / T - Ex+15
where step uses Er+1 2 0 and
/ > 0. This leads to:
ET. E < (1 -1- (V2 V 20 Σt, we
E, + E W/
E, + 2pr,
step @ uses the fact that (1 V V 1 and (1- / v when 1; step 2
uses Incqualities (15 Letting T + 20, we conclude this lemma.
Lemma A.10. Assume that [D.]=+] a(D,-1 D.), where t,a 0, and (Dr) is nonnega- live sequence. We have: Dr < (TT))
19

---

Published as a conference paper at ICLR 2026
Proof. We let K > 1 be any constant. We define h(s) =s---, where > 0.
We consider two cases r' h(D.)/h(D,-)).
Case (1). r S We define h(s) - .- We derive:
IA I
1
a(D;-1 - -D))
a(Dr-1 Di) (1
ax axSD. Dr-1 1())da
ax 1) - k((D))
E ax -.., - (Do-1l-))
where step @ uses [D,] aD D.): step @ uses h(D)) S kh( (hh(D-1)) step @ uses the fact that h(s) is a nonnegative and increasing function that (a (a-b)/((a) - Jn h(s)ds for all a, E [0, oc); step uses the fact that Vh(s) = h(s); step @ uses the definition h(.). This leads to:
DD1 - 1 RD
(17)
Case (2). r¹ >> We have:
h(D:) xh(D.-1)
- [D,]-(+t1) R (D:-1]
$
FTI KTTI
(1+t))-
D > K-T [D,-1
(18)
where step ① uses the definition of h()); step @ uses the fact that if > > 0, then a² > b² for any
exponent + a _ E (0, 1). For any 1 2 1, we derive:
7
[D:] - [D-1]]
(KT+I - -1)- (D-
(KT+I - 1)o
(19)
where step uses Inequality
(18 step @ uses 7 > and Di-1 S Do for all 2 1.
In view of Inequalities 17 and 19 we have:
[D.]
Dr min( (KT+I 1) [Do]
(20)
A
Telescoping Inequality 20 over t from 1 to T, we have:
- [Do] -
2 Të.
This leads to:
Dr D - < O(T-1/T).
B SOLVING THE SUBPROBLEM WHEN =
This section presents a novel Breakpoint Searching Method (BSM) find the global optimal solu-
tion of Problem (3) when &- = 2.
Initially, Problem
boils down the following one-dimensional subproblem:
min 1 + VVP)) + (VZ), s.t. Ve (Vio1 Vref
20

---

Published as a conference paper at ICLR 2026
which can be further rewritten as:
i E arg min } vvec (V) Qvec Vvec)
+ (VV)), s.t. &
- cos(O) sin(0) Fsin(0) cos(9)
)
where Q e R P E 222 and Z R Given h(-) is coordinate-wise separable, we have the following equivalent optimization problem:
min
h (cos(0) sin( + cos + bsin(()
+ccos + (os + esin)
(21)
where 22 Q:
Q13 Q21
Q31.
0.5 Q33) Q23 R2xx1 and y - [s; ul E 2rx1
Z(2,)), p Z(2.)). u (I.)), x - [r:p]
Our key strategy to perform a variable substitution to convert Problem 21 into an equivalent
problem that depends on the variable tan _ t. The substitution is based on the trigonometric identities that cos = =1/ 1=1/v1++ + tan" and sin () + tam( ta (an)))
The following lemma provides characterization of the global optimal solution for Problem
Lemma B.1. We define ((, 8) & ač + bä + cτ + dčs e5 + h (ūx + sy). and w c e. The optimal solution 9 to (1) can be computed as:
sin)())) E arg min F(c, 8), s.t. (c, E {la s1). (c2, s2), (0, 1), -11} [c,s]
where A
1 x
-1
and 82
- V++("
Furthermore, E and
V+
VI+
v VI+
E_ are respectively defined as:
f E argmin, p(t) &
wtit
hi x
(22)
vI+12
1+f
Vi+
i_ E argminn p(t) & -a-bt
w+dt I+12
h xx-1y
(23)
vI+r2
VI+1
Proof. We define w - c- e, and F(2, 3) - aZ + + oz + des + e5² + h(ox + 5y).
With the identity sin (0), Problem 21 can equivalently written as:
i arg min
((co(()) ()) + sin (0)y) + (o) (0) b sin
+ cos (0) + d cos (0) (i + e.
(24)
We first consider the case cos 0 * 0. By Lemma A.8 there are two possible parameterizations for (cos(O), sin())) in Problem 24
Case a). cos
v (+) and sin(0)
tan Then Problem 21
becomes:
V++aan?(0)
04 E arg min a+tmm(0)b
tnn( [+)an-(0) + h +t-tam(0l
(8)
/ (+)
Setting & tan(0), we have the equivalent problem:
f E arg min atbe + widt 1412
+ h( xtyt
VI+
VI+T#
Hence the corresponding optimal trigonometric pair is
cos(f+) V sin(A+)
i+ (++(E,)=
(25)
Case b). cos(
v 1+) and sin(0)
tan(() In this case, Problem 21 reduces Vittancey
Ø_ E arg min tan(0)6
w+tan(0)d
h( taan
Vi-tan(6)2
I+tan(0) V1+tan(0)2
21

---

Published as a conference paper at ICLR 2026
Again letting t tan 0, we obtain
{_ E arg min -0-bt vi+*²
+
w+dt 1442
+ 1h( (1+12)
Thus, the corresponding optimal trigonometric pair is
- cos(0-) VI+(E_)>* sin(0_) vI+((_)²
(26)
Combining 25 and 266 when cos(O) # 0 the optimal solution 0 to (24 satisfies [cos(0), sin
arg minc.. F(c, s), s.t. [c,
{[cos(0+)) sin(0+)]. [cos(O_). sin (8_])- Including the case
cos = 0, that is, [c, s] E {10, 1]. [0, - 1]}, the final selection rule for the optimal pair is
(O), sin(O)] E arg min F(c,s), 6s
s.1. (c, s] E [cco((0)), sin(0+)), [cos(0-)), sin((-))) [0, 10, -0]-
Note that {cos(0), 0), sin(0) uniquely determines 0, and the objective in Problem ( depends only on (cos(0), (i) for some 9. Thus, it is not necessary to explicitly recover the angles 0 and 0_; it suffices to work with their cosine-sine representations.
We describe our BSM to solve Problem (22): our approach can be naturally extended to tackle Problem (23 BSM first identifies all the possible breakpoints/ / critical points θ, and then picks the solution that leads to the lowest value as the optimal solution i.c., E arg min, p(t), s.t. t O.
We assume y; # 0. If this is not true and there exists y; = for some i, then yx. can be removed since it does not affect the minimizer of the problem.
Finding the Breakpoint Set for h(x) A A||x||o
A||x||u is scale-invariant and symmetric with || + (x|lu ||x]|o for all t 0, Problem (2) reduces to the following problem:
Since the function (xx) -
min p(t)
a+Mt w+dt
(27)
V1+1² T437
Given the limiting subdifferential of the Eo norm function can be computed as 0||4/o
t= (see Appendix C.5). we consider the following two cases. (i) We assume
(0), ebe.
(x + ty); = O for some i. Then the solution & can be determined using = There are
2r breakpoints [xi yi y2 Xar for this case. (ii) We now assume (x ty);
for all
Then A|xx + tyllo = 2r becomes constant. Setting the subgradient of p(() to zero yields:
= Vp(t) = [b(1 + () (a + bf)t] +I+ + [d(1 + t²) (w + df)( to, where to = (1 + 2)-3 2 Since to > O, we obtain: d(1 + t2) (w + dt = - (b at) /1 +
Squaring both sides, obtain the following quartic equation: cat + cyt + czt t co
for some suitable C4, C3, C2 Ci and co- Solving this equation analytically using Lodovico Ferrari's method (WikiContributors),
obtain all its real roots (4,12,...,1 [2} with < j 4. There
at most 4 breakpoints for this case. Therefore, Problem (27 contains at most 2r + breakpoints Xax. E1, f2, ..., {
④ 1 32 You
Finding the Breakpoint Set for h(x) & A||x||1
Since the function h(x) A A||x||a is symmetric, Problem 22
reduces the following problem:
f e arg mimm(ip A a+de vi+e
+
widt 14t2
+ A|x+ttyl1 VI-P"
(28)
Setting the subgradient of p(-) to zero yields: 0 0p(t) = tn[d[1 + (2) (w + dt) (6 at)
v1 t2] + X /1 + t². [(sign(x t)) +f) ||x ty]] where to = (1 t2) We consider the following two cases. (i) We assume (x + ty);
0 for some Then the solution
be determined using 7 x.. There are 2r breakpoints
Mar for this case. (ii) We
y,
Yi Y2 yz.
22

---

Published as a conference paper at ICLR 2026
now assume (x + ty); * for all i. We define - {+ 81 X1. + X2 M2c X2= y1 Y1 v
Y2 ya- ya-
R4rx1, and sort z non-descending order. Given i + 2 for all in this case, the domain p(t) can be divided into (4r + 1) non-overlapping intervals: (-00, (Z1, 22), (Zar, +00). In each interval, sign(x x + (y)
be determined. Combining with the fact that 1 > and Xx + ty|l = (o, x ly), the first-order optimality condition reduces to: 0 = [d(1 + +2) (w + dt) + (b - at) v + + y) + (o, + ty)4], which can be simplified as: -xx] = =[[11 (d(1+t))- -+)))((( -w (w+
2t). We square both sides and then solve the quartic equation. We obtain obtain all its real roots { f2, ... {5} with 1 S j S 4.
(at- v t --A- [(o.
Therefore, Problem (88) contains at most 2r + (4r + 1) x 4 breakpoints.
Finding the Breakpoint Set for h(x) - (x)
Since the function h(x) (x) is scale-invariant with h(tx) h(x) forall & 0, Problem reduces to the following problem:
(2)
E arg min p(t) A a+bb
VI+E? I+1 s.t. x + ty 2.
+
(29)
We define I A y > 0} and J _ fdyy < 0}. It is not difficult to verity that fr + ty 2 0} + {- X7 S t,t S
+ (Th max(- X1) S t S min(- x1) _ ub). When 1b > ub, we can directly conclude that the problem has no solution for this case. Now we assume ub 2 lb and define P(t) - min(ub, max (m, 1b)). We omit the bound constraint and set the gradient of p(t) to zero, which yields: 0 vp(() = [b[1 ( (a + bt)d] v1+
Yı
[d(] + (2) (w dt)(2()] tº, ....,,,. } with 1 S S 4 after squaring both sides and solving the quartic equation. Combining with the bound constraints,
where to = (1 +1+ 2)-2 We obtain all its real roots {f1. t2.,,..
conclude that Problem (29) contains at most (4 + 2) breakpoints (P(h), P(12)... P(Gz), lb, ub) with 1 S < 4.
C ADDITIONAL DISCUSSIONS
This section encompasses various discussions, covering topics such as: (i) simple examples for the optimality hierarchy, (ii) computation of the matrix Q, (iii) complexity comparison between OBCD and full gradient methods, (iv) generalization to multiple row updates, and (r) the subdifferential of the cardinality function.
C.1 SIMPLE EXAMPLES FOR THE OPTIMALITY HIERARCHY
To demonstrate the strong optimality of BS2-points and the advantages of the proposed method, we examine the following simple examples of 2 2 optimization problems mentioned in the paper:
min VESI(2,2)
F(V) & IV lf. with
-1).
(30)
F(V) A VESI(2.2)
min
V - B + VII,,,, with B = (1
(31)
Figure [ shows the geometric visualizations of Problems 30 and (31 using the relation ming (V
minvesu((.2)) F(V). The two objective functions exhibit period- icity with period of 2t.. Within the interval [0, 2T),, each of them contains unique BS2-point, while the two respective examples contain and 8 critical points. Therefore, the optimality condition of BS2--oints might be much stronger than that of critical points.
BS2-points vs. Critical Point: Algorithm Instance Study. We briefly review methods that seek critical points of Problem (30) and demonstrate that they may lead suboptimal solutions for Prob- lem (30). As representative example, we consider the well-known feasible method based on the Cayley transform (Wen & Yin 2013). According to Equation (7) in (Wen Yin 2013), the update
rule is:
x QX', Q [L2 + A (I2 _A)],
(32)
where 7 E R. and Á E 2 is a suitable wem-syymttiic matrix. Lemma A.3 shows that the matrix Q is always rotation matrix. Consequently, if X is initialized rotation matrix with det( Q) = 1, all iterates X remain rotation matrices, which in general do not coincide with the optimal solution.
23

---

Published as a conference paper at ICLR 2026
10
.... .(..F)
30
.....(V8')
(VV()((3") ))
-min(F (VV)), F(V;))
Objective
Objective
25
5
20
15
-5
mit F(V) ||V All
minves F(V) |I Bll + 5||
Figure 2: Geometric Visualizations of Two Examples of 2 x 2 Optimization Problems with Orthog-
onality Constraints with A (-1 ) and B (1
C.2 COMPUTING THE MATRIX Q
x as in 8 can be challenging task because it involves the matrix H E IR "rrr However, in practice, H often has some special structure that enables fast matrix computation. For example, H might take diagonal matrix that is equal to Ll for some L 2 0 or has a Kronecker structure where H = H, & H2 for some H₁ RRex and H2 E RR*x., The lemmas provided below demonstrate how to compute the matrix Q-
Computing the matrix Q E x
H1 & H2. have: Q Q1 X Q2. where Q1 = ZH,Z E k and Q2 = H2D E Reck (LZZ)) & L.
Lemma C.I. Assume (8) is used to find Q. (a) If H
(b) 1f = Llnr. we have Q
Proof. Recall that for any matrices A, B. C. of suitable dimensions, have the following
quality: (A & (C & D) (AC) Ø (B)).
If H & H. we derive: Q è & TH = -Z &U.) H. (Z &U.)) =
Z ( (H,Z*) (H,U))] (ZH,7* (U. Q1 Q2-
(b) If H T.n.. have: Q (ZT & U.) H(Z Us) I(ZT & U.) U2)
LZZZ ) Ix-
Lemma C.2. Assume (9) is used find Q. (a) IfH Hi @ H2, have Q | Q2 I,
where Qi and Q2 defined in Lemma C.1 (b) H LI we have Q L 7 L
Proof. (a) Using the results Lemma C.1 [a,, we have:
(Z 2 Up)
Z
X Up) = Q1 8 Q2 1
||Qu llsp ||Qpllsp
(b) Using the results in Claim (b) Lemma C.1 we have:
(ZT&U.) TH &U,) LZZ
L| Z 12 I.
C.3 COMPLEXITY COMPARISON BETWEEN OBCD AND FULL GRADIENT METHODS
We present a computational complexity comparison with full gradient methods using the linear eigenvalue problem: minx F 1 XX, CX), s.t. X = I,, where C E R is given.
We first examine full gradient methods such the Riemannian gradient method (Jiang & Dai 2015 Liu et al. 2016).. The computation of the Riemannian gradient V VM = CX X
requires On=rr time, while the retraction step using SVD, QR, or polar decomposition demands
24

---

Published as a conference paper at ICLR 2026
O(nr²). Consequently, the overall complexity for Riemannian gradient method is N1 On
where N1 is the number of iterations required for convergence.
We now consider the proposed OBCD method where the matrix Q is chosen to be diagonal matrix as in Equality Q9. (i) We adopt an incremental update strategy for computing the Euclidean gradient (X CX, maintaining the relationship Y = CX' for all t. The initialization Y CX
occurs only once. When X is updated via a k-row change, resulting in X = X + (V
N efficiently reconstruct CX by updating Y Y + CUa(V - IJ - in
(nr) time. (i)) Computing the matrix P as shown in (3) involves matrix multiplication between matrices [x*)aa E Rx* and [X*]]
R=** which can be done in Orkk²). (iii) Solving
the subproblem using small-size SVD takes O(k2) time. Thus, the total complexity for OBCD is
N₂ Onrr + rk2 ):3), with N2 denoting the number of OBCD iterations.
C.4 GENERALIZATION TO MULTIPLE Row UPDATES
The proposed OBCD algorithm can be generalized to multiple row updates scheme.
Assume that n is an even number, and & As mentioned in Lemma 2.3 when (9) is used to find Q, the subproblem V arg migi X', in Algorithm reduces to:
min (V. (V h(VCgX)). VESt(2.2)
(33)
One can independently solve (n/2) subproblems, each formulated as follows:
minvesu2.2) (V,
(V, (Vf(X*)
X
( + h(VUg') with B = [1, 2].
minvesup.2)
(V/( (V,
X )s) X with B = [3.44.
minvesu (V, (Vf XX )RR) + AVV with B = [n - [n].
This approach, known as the Jacobi update in the literature, allows for the parallel update of rows of the matrix X.
Notably, one consider & A |B| 2 when h(-) 0, and the associated subproblems can be solved using SVD.
C.5 LIMITING SURDIFFERENTIN OF THE CARDINALITY FUNCTION
We demonstrate how calculate the limiting subdifferential of the cardinality function h(X)
|X o. Given that h(X) = |X] lu is coordinate separable, we focus only on the scalar function
h(x)= lz/o. where ozlo 0. = 0:
}.
else.
The Fréchet subdifferential of the function h(x) = oolo at æ E dom(A) is defined as Ōh(x)) & 2/1 inf2/2
{E IR lim=-z inf
(E,2-s) > 0}. while the limiting subdifferential of h(x)) ---
dom(h) denoted as (x) A {E R Ex + x, h(xx)) - h(x)) E ōh(x*) + E, We con-
sider the following two cases. (i) x + O. We have: dh(x) = (E E R lim inf 2 s-
((,--1) 2
0} = (0). (i) 0. We have: dl(()) 1E E R. limz-ra inf znd =lo 16--6_=-1)
0) (g E R:
lim,- inf 2 1- (5.=)
R.
We therefore conclude that [8] X [|)), E X = } for all i [m] and j [r].
? (0).
25

---

Published as a conference paper at ICLR 2026
D PROOF FOR SECTION 2
D.1 PROOF FOR LEMMA
Proof. Part (a). For any R and (B. we have:
x+ - x
0 U+Us(V
X
UV I+)U U.(V + (V-(V
X (VV I 1-+VV-- -) +)+ (V- -b) U
ø X [(V - -I, V -))++ +V) [+) U
X-Ub(V- - -
-
U-U,(-I: VVVJ,
=
@ U O
0,
where step ① uses X = X UuU - -V step @ uses U3 Up Ig; step uses V I.-
Part (b). Obvious.
D.2 PROOF OF LEMMA 2.2
Proof. We define X á x b,) Q (z & & @ Up). and - U
Part (a). We derive the following results:
X
X|
|V - I
vec(U, - [2)
T Hveo (V - I))
vec(V V (Z) @ Un) J HH @ Uz)ocC( - 11)
IV Is (.ll
V 1x112.
where step ① uses X X
(V - Iz) step @ uses X vec(X) (Vve step 3 vee vec (UUZ) for all R. Z, and U of suitable dimensions; step + uses again; step @ uses the definition of Q.
uses TT R (U)
X vee H
Hve X
Part (b). We derive the following equalities:
x
X
[V - I+) 2
II(V I,)Z
((V T I.), ZZT)
2( V.ZZT) + +(V VT ZZT).
2 V. ZzT +0. +0.
where step 1 uses X X + U. - I Z; step @ uses the fact that V || for any
V E A k step ② uses
(V TV [3) = I - VT_ V + I = 2(I - V) V VT):
step 4 uses the fact that (V. ZZT) = (VT, (ZZT)T) = (VT ZzT) which holds true as the matrix
ZZ is symmetric.
26

---

Published as a conference paper at ICLR 2026
Part (c). We have:
X X
[V X
1A8 I
||Uel II(V - - |
| sl
II - I . 10 |X
1
II
L
1|
2(I V.
where step ① and step @ uses the norm inequality that | |F I ||A|| || for any A and
X; step @ uses | = | |X|| sp for any X St(n,r): r); step @ uses the following
equalities for any V St(kk, k)::
IV - Ix |l
V i + 4 2 V) | 12 L 2 V) 2 V).
D.3 PROOF OF LEMMA 23
Proof. We define K(V; X B) 1/l - L I + (VZ) + (V, (V (V))
[Ra) + ë. where
Z A U. and ē h(U X f(x) (Is, (X)(X((XX) (Ra) is constant.
Part (a). Using the definition of K(VV X B), we have the following equalities for all E St(kk
X(V;X',,) B)
A
8+ V I Qul + (V, (V hoo) +((Z)
8+ V l 1 I || (V, | f ( (BO) + h(VZ)
@ a VI a/ wat (aae(Qvec(eaa QQeec (1)) ))) + 4 l + [1. L - V) + (V, [VX(X' + 4(VZ)
2 a+ 2
V Q + (V, VS(X') D8 mat Q vec (1)) alg) + h VZ) + +
where step uses Lemma 2.2 c) that: |/l I (I. 1 V); step uses the definition of P.
Part (b). We consider the case that Q is chosen to be diagonal matrix that Q = <I,, where defined in Equation (9). Using V E St(k,)), the term V simplifies to a constant with
!! ! = sk. We can deduce from 3
V arg
min VESu(K,
P(V) & (V.P) + ((X).
(34)
In particular, when (XX) = 0, Problem becomes the nearest orthogonality matrix problem and can be solved analytically, yielding closed- -losed-form solution that:
V E arg min (ESt(K,A)
1 V - (-P) PM = -PM(P)=-UV)
Here, ÜDiag(s) VT the singular value decomposition of P with U, V e St(kk k), s R* and 0.
Notably, the multiplier for the orthogonality constraint VTV = I can be computed as: A Vp & (UDiag(s)VTTT [-ŪVT] VDiag(s) ÜUVT & (s) VT 0, where step
① P ÜDiag(s) TT and V UVT, step @ uses Ü Is; step 3 0.
D.4 PROOF OF LEMMA 2.5
Proof. Any matrix takes the form V (ab). The orthogonality constraint implies that
V E (2, 2) meets the following three equations: a² 12, e² + d2, bd.
27

---

Published as a conference paper at ICLR 2026
Without loss of generality, let c sin(f) and d cos with 0 E R. Then obtain either (i) a = cos , b = sin(0) or (ii) a = - cos(0), = sin(f). Therefore, we have the following Givens rotation matrix a and Jacobi reflection matrix Vgrt
Vool
cos(() sin(0)
sin(0)
Viref A
cos 0) sin (0) cos (0s(())
sin
cos(O)
Note that for any b, c, d E R, we have: det( d) = ad - bc. Therefore, we obtain: det(Vg (V = 2 (0) = 1 and dct(VV cos cos2(0) s0) = -1 for any 0 IR
E PROOF FOR SECTION 3
function [Q Jacobi Gi vensQR (X)
size (x, 1); Q en)(n) R
2
for j-1:n
for i=n:-1: (i+1)
4
[i-l;il; V Givens (R(1-1, R(i,, 1));
5
R(B RB, Q(;,B) =);;BB +V;
6
(i--j+1
R(),])<0)
7
V [-1 0; C -il; $ or = [-1 0; 0 11:
8
:) = -BB, ():,B) = -); VV;
9
end
10
end
11
end
12
if (R k)
13
V = [1 0,0 -11; B(B = VBBB(, Q(:,B) = -(;,B)+V;
14
end
15
16
function Givens (a,b) $ Find a Givens rotation that v' [aab]] - [r;0] if (b--0)
17
18
19
20
else
if (abs (b) > abs (a))
22
tau = -a/b; 8 1/sqrt (1+tau"2); = s*tau;
23
clse
24
tau -b/a; с 1/sqrt (1+tau"2); s c*tau;
25
end
26
end
27
v
[c 8:-8 c];
28
Listing Matlab implementation for our Jacobi-Givens-QR algorithm.
E.1 PROOF OF THEOREM 3.1
Proof. Part (a). First, recall the classical Givens-QR algorithm, which is detailed in Section 5.2.5 of (Golub & Van Loan 2) This algorithm can decompose any matrix X E 2ax= (not neces- sarily orthogonal) into the form X = QR, where Q is an orthogonal matrix (Q E St n)) and R is a lower triangular matrix with R,j 0 for all i < j. achieved through C m(n-1) Givens
rotation steps.
Combining the result from Lemma A.5 decompose any orthogonal matrix into the form = QR, where Q E St(n, n)) and R diagonal matrix with R. (-1,+1) +1} for all E [n].
conclude that classical Givens-QR algorithm can
We introduce modification to the Givens-QR algorithm, resulting in our Jacobi-Givens-OR al- gorithm presented in Listing [ This algorithm can decompose any matrix X E St(n, n)) into the form X = QR, where Q and R I,, using a sequence of C Givens rotation or Jacobi reflection steps.
28

---

Published as a conference paper at ICLR 2026
Please take note of the following four important points in Listing
a) When we remove Lines 7-10 and Lines 13-15 from Listing it essentially reverts to the clas- sical Givens-QR algorithm. Givens-QR operates by selecting an appropriate Givens rotation sim(0)
matrix V = [_s((()) cos(0) with suitable rotation angle to zero-out the matrix element n)) and bottom top (i
Rij systematically from left to right (j
1)) within
every pair of neighboring columns.
b) Lines 7-10 and Lines 13-15 can be viewed as correction steps to ensure that the entries Rj for all = n1.
c) Line 7-10 is executed for (n - 2) times. In Line 7-10, when Javeni-G detects a negative entry 1,-1, with i j 1, it applies a rotation matrix V _ ( - to the two rows B = i - 1. to ensure tha 11;-1,11111
d) Line 13-15 is executed only once when det(X) -. In such cases, we have R (D - and (Ræe) = -1, 1, where
is the two indices for the final rotation or reflection step. To ensure that the resulting Ras is an identify matrix, Jacohi-Givens-OR employs reflection matrix V = ( i ), leading to V 12-
[n-1,, n,
[
Therefore, we establish the conclusion that any orthogonal matrix X E St(n, n) can be expressed as D We ... W.WWW where W; Us, V.U + Ugg U... and Vi E St, 2) is suitable matrix associated with B, Furthermore, Ví, V; = I2, we have Vi, W; = I,, leading to D = In. This concludes the proof of the first part of this theorem.
Part (b). For any given X E St(n, r) and X E +St(n, r), we let:
m)) (x(x")),,
(35)
where Pst(n.n) (Y) denotes the nearest orthogonality matrix to the given matrix Y.
Assume that the matrix X[X has the following singular value decomposition:
X(X)) = UDiag(z)VT, E {0,1)", U E St(n,,), V E St(n,n).
Therefore, we have the following equalities:
Diag(z) U'X(XQTV.
(36)
D UV St(nnn).
(37)
Furthermore, we derive the following results:
z E {0, 1)
Diag(z)T = Diag Diag(z)
UDiag(2)) Ding(a)Ding(a))1UXX
UVTx Xx
1 -
where step ① uses (36); step @ uses UUT I,, VVT = In, X I,,, and Xx I,: step uses 37 We conclude that, for any given X St(n,r) and X St r). we can always find matrix D St(n, n) such that o X.
Since the matrix D E St(n, n) can represented as D = WWWW where W. = UB. V.UT By for some suitable V; E St(2, 2) (as established in the first part of this theorem), we can conclude that any matrix X St(n,,) can be expressed as DX WW,
'Alternatively one can use the reflection matrix V (-10) instead of the rotation matrix V & that R. R 1.
0
29

---

Published as a conference paper at ICLR 2026
E.2 PROOF OF COROLLARY 3.2
Proof. We denote e as the i-th canonical basis vector in IR".
We denote the set (Bi, Bc, } all possible combinations of the index vectors choosing
items from n without repetition.
Part (a). Fix any 2. By Theorem 3.1 for the case 2, for every D St(n, n) there exist
index pairs (pi, 9:) and matrices vy!a) E St, 2) such that
D = w ... W (
where
W!!i I, + U (2)
(250)
U lep Eg.l E IR
We let
V; _ (222) 0
kt) w I, UB. (V: I U
Tx-2
By construction, Wj acts as 2/22) on the two coordinates Pj, Oj and as the identity on all other coordinates, hence W3 = w as linear operators on R Therefore
(2) W Weg W. D W
which proves the first part of this corollary for any k 2 2.
Part (b). similar argument to that used in the proof of Theorem bb] yields the second part of this corollary.
E.3 PROOF FOR THEOREM 3.8
Proof. We X, X, and X to denote global optimal point, BSa-point, and critical point of Problem ( respectively.
Setting the Riemannian subgradient of K(V; B) w.E.t V to zero, we have 0 E @m/ X, B)
G(V) θ VVG(v) where G(V) = a(V I.) + U mmat(Hvec(X* - X)) + Vf(X) +
0h(X*) U and X X + (v(V - I.) Letting V Ic, we have the following necessary but not sufficient condition for any -Sx-point:
VB E ( 0 U- (GXT XG")Ua, with G GEVS(X)+ + (X))).
(38)
Part (a). We now show that (critical points X} 2 (BSk-points X} for all > 2. We let G
v (X) (k) Using Lemma A.l have:
On GXT - TGT
(0 X) = - -
Op, = G G
(39)
X 0n, = XT X
$ 0. XTG - -
On,n X(XTG - cTXRT
O X x
GGT
AG
where steps C and @ that the necessary condition in Equation 38) equivalent to the optimality condition of critical points.
X I, step 3 uses Equality (39) that G X We conclude
30

---

Published as a conference paper at ICLR 2026
Part (b). We now show that -Bointss X} iglobal points X} for all k
We define (V) - X Us(V L and B)) _ f(X) (V L. ( as) + 2 Q
(2.3.... ..... n..
Q+al. + +h(Ug. + XVV. We let V( E St k) and
B;
(Biter 1=
C
We derive:
K((a; X, B:), VB;
F(X) h(X) f(X)
4(X) + f(X), VX rt(nrr)
4(X + (V -4JUB + f(X + UB, 1. X), (V AB,
M(X,, (V + ((()),(V()) W((),, VB,
K(V (0) X. Bi), VV VB;
1 minvesu(k.)) KC X, Bi), AB,
(40) +
where step uses the definition of (CVV X, B) & | a + +A(UT_X) + 3 X); step @ uses the definition of X; step @ uses the basis repre-
-I, (xa))
sentation of orthogonal matrices for all 2 2, as shown in Corollary 3.2 step @ uses the definition
of (V); step & uses the same strategy in deriving Inequality (10) This leads to:
I E arg min (k,)(
K(V;; B,),) VB,
The inclusion above implies that (BS&-points X 2 {global optimal points X).
Part (c). We now show that (BSx-points X} 2 {BS -pontss X It evident that the subproblem of finding BSt-pooint is encompassed within that of finding 155t+1-points stationary point. Thus, conclude that the optimality of the latter is stronger.
Part (d). The inclusion fcritical points X}- C (BSg-points X} may not always hold true. This can be illustrated through simple examples of x 2 optimization problems under orthogonality constraints (see Appendix Section [.T more details). Lastly, it is also evident that the inclusions (BS2-points X} S oglobaloptimal points X) and (BSx-points X} not always hold true.
(BSk+1-poing X) may
F PROOF FOR SECTION 4
F.1 PROOF FOR THEOREM
Proof. We define K(V; X B) A / Is I1 +h(VZ) (V, (V, (XXX + ë, where
Z A X and ë h(U Xt) f(X) (I,, VXXX') is constant.
We define ē 2 (F(X") F(X*))
Part (a). First, we have the following equalities:
h(x+1)) h(x') _
+ Up - Up- U.
- M(U V + M(U,- U. X A(U,U - h U U
m MVV' x') x(U
(41)
where step uses X U. V + Uge UJ X' as in 4 and I = UgU + Ug- UT; step 2 and step @ use the coordinate- wise separable structure of h(-).
Second, since V E arg minvee((a.)) X(V; it follows that KC(V: x B) < ICI This further leads
h IVV. x') + 1 V' - -Ill@aale (v' - I (X+T (B)) S ((UX'')
(42)
31

---

Published as a conference paper at ICLR 2026
Third, we denote X X3(V") and derive:
-(x+t)-- -
(X)(V9 x (V)) 1
(U,(V' - I+)
2 V
1
(V' - L I,
V
130 + V -4
(43)
where step @ uses Inequality (2
step @ uses Lemma 2.2 a); step 3 uses Q Q
Adding 41 and together, we obtain the following sufficient decrease condition:
F(X(1)) F(X') - 2 V I -
-
(44)
where step 1 uses Lemma 2.2 c).
Part (b). We assume that B² is selected from {B, j randomly and uniformly.
Taking the expectation for Inequality by each iteration:
we obtain a lower bound on the expected progress made
Eg [x F(X+) < E [ V IIll]. Telescoping the inequality above over 0, 1, ...., T, have:
Egrli E IV - [+]] < Ber (F(X)) - F(X <
F(X
where X denotes the limit point of Algorithm such that
As a result, there exists an index & with S
Egr [I[ LI a(7+1) -(X*)) T+1
F(X*)]
(45)
Furthermore, for any t, V is the optimal solution of the following minimization problem X V E arg miny KC(V;* B²). Since V is random output matrix that depends on the observed realization of the random variable B', directly obtain the following equality:
dist(LR, arg miny XCV;; X B,))² Ee: V - -.
(46)
Combining (45) and (46). we conclude that there exists an index with [0. 7] such that the associated solution X' qualifies as an e-BSt-point of Problem (1), provided that is sufficiently large such that 2
T+I S c.
F.2 PROOF OF LEMMA
Proof. We define & e B as the element-wise subtraction between sets A and B.
We let H++ E (
X 1+),, and define:
So U.. V (xx H X U RAXE
(47)
Ω (
+
T Uge E g
(48)
S2
A (V(')) V][X+ e g
(49)
Part (a). First, using the optimality of V for the subproblem, have:
Og = G T
where G ma + aI vve -1)) + + U VX(X + H H "X')"
ATi
Using the relation that G = Tı + Tz, we obtain the following results from the above equality:
0g (T, + T2) - v' + Tv
02 Yi + 24 22 - V' + +22)
Q v'(,, + + + - - -2
(50)
32

---

Published as a conference paper at ICLR 2026
where step T uses Y₂ 92 + f22-
Second, since both B' and B are randomly and dependently selected from {B,) with replace-
ment, each with an equal probability of
ج for any A E Ruxn we have:
Bg [| g.+ AU | ch E AUB.IE Eg U AUg
Using the definition ç² _ (B ,B.,...,B),, we have:
Ee U 1l Eee B
(51)
Third, we derive the following results:
Bee+n [dist 2 (0, A X =
Ege+ +
e X 11 Ugeen
Re (0F() IJ
x x ()
IN Eg [1l20
12
8Eq (119
P li + 2 [| 12
1
SEg
[1
1 + g (T, + .2 (2)
22) T -T, f2 11
8Egg (1
2 + 6. V V 1 6E.. V 1
(IV
22 l2
(52)
where step ① uses the definition of a X 3(+1) the point V I step @ uses Equality
step @ uses the definition of Qg in Equation 47 step @ uses Lemma A.2 and the fact that (a + b) < 2n + 2b for all a, b E RR step
51 with À )F(X++1(XX+1)
X (aF(X) (
uses Equality (50): step @ uses the inequality (a + + c) S 3a² + 362 + 3c² for all a, b, c IR.
We now establish individual bounds for each term in Inequality 52 For the first term 8Eg [| 2 in (52).
we have:
8Eqe Og ( 8E. [| (V H x -X
[ V
H I2) x U.
8C7 V I 11
(53)
where step @ uses X X (v I
[ step @ uses the inequality |XY|| s
| so for all X and Y repeatedly, and the fact that | S Cp for all X E r)( and all G GE OF())
For the second term GEg
V'rIV'
- Y | in 52
we have::
6E.z V V Till
12
V
+ 12E (1
24E [|
8
24 |||mat (iQat((Q-aa + aIg 11)
r
24 [|| + a V| - [+]]]]]
^
24(L, + +a² Eg V
1)ll
(54)
where step uses the triangle inequality; step @ uses the inequality |XY|| < |X|| Y| for all
X and Y; step @ uses the definition of Ω2 in (48); step @ uses the fact that [Q] sp T.g.
33

---

Published as a conference paper at ICLR 2026
For the third term 6Eçz V V - 2 in 52 we have:
c [ 2 v 22 [i]
Wn 11
6 [ (V' - I+)) + (V- 10]
12B [||1] 1pp | Is [|
12E
V
sp V I
E.
V
(55)
where step ① uses the fact that V + V'Ω = 0; step uses the norm incquality; step I uses the fact that ||21ls || fx [V + I|F which can be derived using the norm inequality.
+ ls 2 [IV + IIssp I
For the fourth term 6E (1) V V 221l2
52
we have:
6R
1211
I^ I^o
12E. V V | + 11
24 (1
||
24E [1 [V
I
24E (IV(()) - Vf(*+*)
I^ I/8
R4 |||X -x--11
24L
IIV''
V Is 112.
(56)
where step 1 uses the triangle inequality; step @ uses the norm inequality; step @ uses the definition of 22 = UV/ V/[X+>>||X]] in 449): step 4 uses the norm inequality; step & uses the fact that (X) is Liipschitz continuous; step @ uses Lemma 2.2 c).
In view of 53
54
55
56 and
52p we have:
E Ee++
110M
K. (I; X + || $ (a
+
02 + + c4) Be (' V
Islle),
where 0 = 82 2 c= - and c4 = 241?. Defining Φ 2
we conclude that | [11 B/+1, || v &c.: [ - I.
Part (b). We show that Eg /dist 2 (0, (x))( 7 Eg [dist (0. (I B'))], where
7 2 0// -O For all D' OF ( x10F(X')) we obtain:
= E
2 (D) 1 12 +
@
E,E,
2
| B:
D'U l|
*
E [ D'U
E.
||| D Ugr [|
(57)
where step @ uses the fact that D = 0 for all E [n]: step @ uses Claim (a) of this lemma with
D = 0 for all [n]: step @ uses Eg [|UT
WUg ||
& - |U] WUR ll as B' are
34

---

Published as a conference paper at ICLR 2026
chosen from {B. randomly and uniformly; ① uses the definition of y. We further derive:
x
yB (OF( | e x'
(58)
where step @ uses the definition of (f )); step @ uses X x
x Isi step @ uses the inequal-
ity that || 2 < | for all X E St(n, r); step 4 uses Equality 57 step @ uses the definition of OMK X )").
F.3 PROOF OF THEOREM 4.6
Proof. We derive the following results:
Eer [distt (0, 0..F0FF22 (0
We i^8
7 Egr+ |dist (0, [IR X T+1, +T)]
7. Φ Eer ||V 11
7.06 741
where step uses Lemma 4.4 b); step @ uses Lemma 4.4 a); step @ uses Inequality (45).
Therefore, conclude that there exists index i with ē [0, T) such that the associated solu- tion X qualifies e-critical point of Problem (1) satisfying lEee (dist" F0 F( S
è provided that T sufficiently large to ensure 7 Ø- T+1 S.
F.4 PROOF OF THEOREM 4.10
Proof. By Theorem 4.2 a) and Theorem 4.6 the composite function (X & F(X) + (XM is
monotonically non-increasing, i.e., F, (X s FX(()). Moreover, the sequence (X has limit point X
Since (X & F(X) + (X) a KL function by assumption, Proposition 4.9 implies that there exists an index t, E N such that, for all t 2 t.,
(' F( & dist(0,OF,(x)))
(59)
Since vp.)) is a concave desingularization function, we have: yo(b) + (a (a--- (a) S s(a). Applying the inequality above with a = F(X') - F(X and b = ( F(X we have:
(F(X)) F(X++
F((x9))
s ((F(X') - F(X - - F(X*)).
(60)
35

---

Published as a conference paper at ICLR 2026
Part (a). We derive the following inequalities:
(B++1)² A Eg. [I] [x]l
Re (e((F(X)) - [(X
Rei [(pr - ((+1) ((F(X)
I^@
2. Ree - 9(41) dist(0, (F,()'))
ele I^ Ie I
2 1g::((e 9441) 110mF(X))ll
2. ((0 9411 X B'IlE]]
2 Egr [(p: ( 10 Lllel
Vo 70 (P++1) Et:
where step S uses the sufficient decrease condition as shown in Theorem (.2 step @ uses Inequality 60); step @ uses Inequality (59); step @ uses Lemma A.7 step @ uses Inequality (58); step @ uses Lemma 4.4; step ? uses the definitions of {MC, (Pc, E.).
Part (b). Applying Lemma A.9 with Pt Kipi with Pt 2 Pe+1, for all i 2 1, we have:
j= Ej+1 S Eg + 2p.-
Using the definition of D, & ΣØ= Ef+1 and letting t. we obtain:
D. S E + 2p:- ① Eg + 2kpp I Eg + 2kpp S 2/ + 2kP1.
*
where step ① uses Pt Kip; step X uses f S P1: step @ uses E, (Eg [|| - I [IF)) 1/2
and | IIllF - | upper-bounded.
F + I S V k Vk. We conclude that D, Ejer Ej+1 is always
Using the fact that |X - X <
V 1 as shown in Lemma 2.2 c). we conclude that
Σ Eg Xi+1
-X
X []] is also always upper- bounded.
F.5 PROOF OF THEOREM 4.11
Proof. We define sp: _ *(s*), where s A F(X') F(X*)).
We define E,+1 - (Ee( V - Ibll)) 1/2 and D; = Ej Ej+1.
We have: D = E, 2F.
First, we have:
9
X X
i/ m
Ejr X Xi-|||
J= V I F
Ejer Ej+11
Dr:
where step ① uses the triangle inequality; step 2 uses X X V Lemma 2.2 step 3 uses the definition of E.+1: step @ uses the definition of Dr. Therefore,
ILl... shown in
suffices to establish the convergence rate of Dr.
36

---

Published as a conference paper at ICLR 2026
Second, obtain the following results:
E
W/@ 1/8
F(X))I)
B X B')I
IV 1
70 V ++lF
Vy
661)
where step ① uses uses Proposition 4.9 that d0st(O, OF, (F
step 2
uses Lemma A.7 step @ uses Inequality
(58): step ① uses the Riemannian subgradient lower bound
for the iterates gap in Lemma 4.4, step @ uses the definition of E, 2 [/ - I [
Third, using the definition of Dr, derive:
Dr 1
E i=t Es
Eg + 2x4pt
E, + 2xc [[s"]r]] 1-a
1-a
E, + 2xc fc(11 -0) (p()))
E, + 2xc {e(1 -a) V 1-a
D:-1 D, 2xc {c(1 σ) /7 D2-1- D.)} 1-a
Dz-1 - Dr 2xe [c(1 - (o(1-o)v Tc 1-0
1-a
=
(D=-1 - -D+}
(62)
AR
where step ① uses Eis Es41 S E, + 2kpL, shown in Theorem 4.10 b); step @ uses the def- initions that Pp & ((s"), and (ss) cs--a; step 3 uses (s) c(1 σ) [s]_ -0 leading to [s*]" = c(1 - o). (st) step @ uses Inequality 611; step @ uses the fact that E, = D D..
We consider three cases for (0,1).
Part (a). We consider o = 0. We have from Inequality 61
0 < Eel [( Vo
Ree I- c(1- a((-a))169-> v Vo
Eee [s]-"; step @ uses 0 and E, = D Dr.
(63)
where step ① uses ('(s) c(1 0)
Since E, + 0, and 7, o, > 0, Inequality 663 results contradiction E, 2 cv/70 > 0. Therefore,
there exists t' such that D, = for all t', ensuring that the algorithm terminates in a finite number steps.
Part (b). We consider a E (0, 1) We define A 1-a 2 We have from Inequality 62
D < D.-1 - D, + (D:-1 D.)
D,-1 Dt + (D,-1 D:) -
Di-1 Fu 1,x42
75-1.2+1
(64)
where step D uses the fact that x < x . -l for all a E (0, 4., and x = D D. E [0,7]. Therefore, we have:
Dr < Di
x-1 R+d
T-1
50-1.+2
37

---

Published as a conference paper at ICLR 2026
Part (c). We consider a E (2 1). We define w 1-a (0, 1), and 7 - 1/w 1 (0) We
have from Inequality 662
Dr-1 Dr +R. (D,-1 D₂) L-0
De s
_
R(Di-1 - D1) + (Di11 - D.) ". (E,)
S x(D;-1 - D + (Dr11 D.) -
(Di-1 DI)
t. (R (
AR
where step uses the definition of w and the fact that Dz-1 Dj Egg step uses the fact that (0.. x--w S -l- w if w E (0, 1). We further obtain:
D D.j 1/ S (D2-1 D2) k//w
=[Dij+41
Applying Lemma A.10 with = k//W,, we have:
Dr < O(T 1/7) @ O(T 1/w-1) O(T 1-o ) OTT 20-1)
where step uses A 1/w 1; step @ uses w A 1_g.
G ADDITIONAL EXPERIMENT DETAILS AND RESULTS
This section provides additional experimental details and results for our proposed methods. We first introduce nonnegative PCA as an additional application, describe the datasets and experimental settings, and specify the compared baselines for é f1-regularized SPCA and nonnegative PCA. We then report extended results on fo-regularized SPCA, f--regularized SPCA, and nonnegative PCA, demonstrating the effectiveness and robustness of our algorithms across these settings.
G.1 ADDITIONAL APPLICATION: NONNEGATIVE PCA
Nonnegative PCA is an extension of PCA that imposes nonnegativity constraints on the principal vector (Zass Shashua 2006 Qian al. 2021). This constraint leads nonnegative represen- tation of loading vectors and it helps to capture data locality in feature selection. Nonnegative PCA can formulated as: minx (CX, X), X where C R"x* the covariance matrix of the data.
G.2 DATA SETS
To generate the data matrix A E mx we consider 10 publicly available real-world or randomly generated data sets: 'wla', 'TDT2', '20News', 'sector', 'E2006', 'MNIST', 'Gisette', 'Caltech' 'Cifar', 'randn'. We randomly select a subset of examples from the original data set. The size of A mmxn is chosen from the following (m27)) 300), (500, 1000), (8000, 1000), (6412, 1000), (2000, 1000), (60000, 784), (3000, 1000), (1000, 1000), (500, 1000]]. We scale the matrix A to have unit Frobenius norm by setting A = ||]r and let = ATA E R
A
G.3 ADDITIONAL EXPERIMENT SETTINGS
Compared Methods on L Regullarzd SPCA. We benchmark OBCD against the following
state-of-thh--ar algorithms: (i) Randomized Submanifold Subgradient Method (RSSM) (Cheung
et al. 2024); (ii) Linearized Alternating Direction Method of Multiplier (LADMM) (He & Yuan.
2012: (iii) Riemannian Subgradient Method (RSubGrad) (Li et al.] 202T]; (iv) ADMM (Lai &
Osher 2014 (v) Manifold Proximal Gradient Method (ManPG) (Chen et al. 2020). For RSSM
38

---

Published as a conference paper at ICLR 2026
and RSubGrad, the subgradient G' aF( iterate X is taken G' CX* + Asign(X'), since sign(X) is a valid subgradient of ||| All competing methods are initialized with a random matrix, producing five variants: RSSM (mnd), LADMM (md), RSubGrad (md), ADMM (rnd), and ManPG (rnd). For OBCD, we employ random working-set rule with identity initialization, denoted by OBCD-R(id).
Compared Methods on Nonnegative PCA. For Nonnegative PCA, we compare OBCD with
two leading infeasible approaches: (i) Linearized ADMM (LADMM) (He & Yuan 2012 Lai &
Osher 2014), (üi) Penalty-based Splitting Method (PSM) (Yuan 2024 Chen 2012), and (it) Rie-
mannian ADMM (RADMM) (Li al., 2024a). Since LADMM, PSM. and RADMM are infeasible
methods and may violate the nonnegativity constraints, we evaluate the quality of intermediate solu-
tions using a surrogate objective, ((X) + 1000]| min X) ||F with X St(n, r), which penalizes
any violation of feasibility.
G.4 ADITIIONAL EXPERIMENT RESULTS
Results on Regglaular SPCA. For À E {10, 50, 100, 500), Figures present the con- vergence curves of the compared methods In-regularized SPCA. Across all setting, OBCD-R consistently achieves lower objective values than competing methods, further reinforcing the con- clusions drawn in the main paper.
Results on Le-Regularzedd SPCA. For À E {10, 50, 100, 500). Table and Figures 7 10 report objective values obtained by all methods with
20. Two observations follow. (i) ManPG is generally faster than LADMM, ADMM and RSubGrad, which aligns with the findings reported (Chen et al. 2020). (ii) OBCD-R consistently achieves lower objective values compared with (LADMM, ADMM, RSubGrad, ManPG), demonstrating its superior solution quality.
Results on Nonnegative PCA. For {10, 20, 40, 80), Table 3 reports objective values and feasibility violations measured by min X)IFF. while Figures 11 shsoww
show the surrogate objective f(X) + 1000|| min X)IIF. Two key conclusions can be drawn. (i) The proposed methods gener- ally achieve the best overall performance, and OBCD-R often substantially outperforms LADMM, PSM, and RADMM by locating stronger stationary points.
39

---

Published as a conference paper at ICLR 2026
PALAARINH RADARRAIA
Relaiiw Oojecteee
RADARAAGA
Pelastee Orgendres
Peblltv Orgectie 8
Relatee Orgerdaee
15
Time (seconda)
Time (soconds)
Time (seconds)
Time (S0conds)
'wia""
(b)
10
ARMBRRT POOEISE
ARAMINAAI PICVL
Theaaive Ocpecios
Daetaias Obpectiv
Thladive Ohvectvv
Detaiaw Obpectivs
PISELO
100
10
60
(c) '2000
'MNIST
'Gisette"
(h) 'CunCullech'
Thbaiive Chjecivve
BADABALL)
Pebattee Ohierton
BA BADABBOO
10t
10"
ORCD-RUA
10"
Time (secorde))
Time (secorde)
'Cidae'
① 'randa'
Figure 3: The convergence curve for solving Lo-regularized SPCA with À = 10.
40

---

Published as a conference paper at ICLR 2026
1gP
ARMMIRT POOESSS
Deearive Okjectos
Duniiis Olypctio
PREMUCT
Debiii Olpectov
PIEVE
Relriww Obyective
PEMUO
RATIRRAAN
DBCD-301
(a) 'wla'
(b) TDT2'
(c) '20Newss'
'sector'
10*
Reaatr Olyasion
RADADNBLA
Reeater Clyysioo
RADAWALA
Reatee Chiesion
LADAINOOD
DOCD-RiA
TM"
Realee
OBCTD-900
Time (secorde)
Time (seconde)
Time (secorde)
Time (secordd))
'20006'
"MNIST"
'Giseme"
(h) 'CunCaltech'
1nt
10t
Releste Okjjetivs
Peleto Osjecivs
3
TRICT-RIAI
10"
Time
(soconds)
Time (Soconds)
(i) 'Cifee'
'rando'
Figure 4: The convergence curve for solving Lo regularized SPCA with À = 50.
41

---

Published as a conference paper at ICLR 2026
LLASMOOA RADARNAGA
Relaiiw Oojectire
Palatre Orgendan
RADAROOIA
PAAAlt Orgectve 00
RACIMNAG4)
RACIMNAS
-0
5.
Time (seconda)
Time (Soconds)
Time (seconds)
Time (soconds)
'wia""
(b)
"soctor"
10
ARMBRRT POWEIS
ARRMINIAAI PICVL
Thlaiiw Ocpecios
Deeaain Obpectiv 15
Theadive Oloectov
Netetww Ofpectiv
PEMLIO
RATIORAANE PIVAINO
R. RATRRMANO
&.
10
da)
(c) '2000
(f) 'MNIST
'Gisette"
(h) 'CnnCullech"
Thbaiibe Oljecivv
BADABALL)
Redattve Oljerton
RADABAD4)
10"
(secordee)
(secordd))
'Cidae'
'randa'
Figure 5: The convergence curve for solving L Leguulaarr SPCA with À 100.
Nweate Obsecivv
Thratres DIUNSTAN
PROAGC
10
LASMMOUd
DOCD-2RAO
CHCCD-RR
(secorde)
(seconte)
(secordd)
(a) 'wis'
(b) TD12'
'20Ness'
(d) 'sector'
Pelaiiv Oojectin
Pelalive Ospeches
Peletie Ocjecten
RATIMMAGES
Peleiiee Oqgectes
EL DATRWAAGEE PISEUO
RADARAREE
RECD-RIS
DRCD-RiA
10"
10
20
(seconds)
30
60
(Soconds)
30
43
50
40
10
22
30
40 (Seconds)
30
'E2006'
'MNIST"
'Giseme""
(h) 'CinCaltech'
Thetrivw Olopectes
PIMUUC
Objecive
10¹
70"
LECD-BDAl
10*
107
10
Time (seconds)
Tima (secords)
(i) 'Cifa'
ii 'randn'
Figure 6: The convergence curve for solving Legreuularized SPCA with À = 500.
42

---

Published as a conference paper at ICLR 2026
BSSM
LAUMM
ADMM
MankG
LADMM
Manto
211
30 200
1016 50
127 0011
393375
645.55m
159
WIP 700
T8R0 191
0017 0379
1078 711
1503 105
TDTT -00 1000
4796 905
997
376,055
2796 315
199
159636
TDTT
24811 390
1142 577
5589 161
130M 188
599
959/663
1E290 6412- 1000
5088 699
211.55 538
252.331
2548 919
199.900
199:545
1E580 020mm-6132 12 1030
25683 861
1005244
9256-712
13085 491
199 0939
999:334
EAXe 2000 (000)
4201 1944
201 023
240.955
2873 292
203.000
199.341
L2K6-2200 1000
29665 881
1085-355
4202-980
15000:413
1000.000
9933333
MNISS 50000-784
334,146
3077.66
3335.519
090-992
CraCidtech 3000- 1000
5274 750
203.177
237.583
2962 145
199.990
193:331
CxsCSdtech3000 1002
26-42.906
1075 922
5548.585
14435 079
(50) 0900990
099:335
Cila- 1000- 1003
5526 600
199.900
451.880
3007 105
190.990
199:33
Cile- 1000 1000
28174.155
1101-222
6080.340
14823 673
090.000
5903285
wade-500-1000
207.00
257.307
2909.599
rands 5200
25917.437
390
4536.156
581
199.000
dra-mo
RS5M
modi
LADMM
(rnd)
ADMM
MasPG
[nxl)
OBCD-R
(40
dep-m-n
RS5M
LADMM
(ma)
(ind)
ADMD
(md)
MaaPG
OBCD-R
(1)
26,
immint
20.
Trre immint
TED1 500, 1000
46043 565
2210:215
1:221 6-2
1992.993
1999/36
T1012 2D1128003000
2000
245899330 363
11036.53 202
173800 198
15/885.090
9996626
2choow 2agg- 10000
53122 247
2204
22141.588
10061 607
1993997
1999.513
5903 1000
2575132913
101.9.354
123658.11 121
1523413123
5002.00
0989055
secken 541--100
51434
2222
17521.156
623
1992.990
1999.834
14220-5412 1030
9999.915
199887
00000000
07001224
MINSS--5500 650000000000 184
461911009 029
3057 975
13300 305
30005 421
1999.992
1999.940
MUNIST-65000-77
784
220035 432
10306.371
148877-32 726
145566.75
794
0900.0021
Glastic 3:00. 1000
51596 303
2201 300
15971.871
30008: 736
1925.990
1995.989
Giacte 9000-1000
281983 001
10313.000
202913.550
152724405
0902.00
C7000009
9937.896
1992.990
CasCliteds- 3000-1000
45
10415.351
158856.813
09396927
2292.400
Cila: 1000-
362258
caale 500 1000
52275
2309 568
1489).188
30522 549
1995.990
1996.988
saale 500 1300
257825 41
10229
80831,264
1570500 52.
1993.999
9999.955
Table 2: Comparisons of objective values for L1 regularized SPCA. The 1st 2,2dd and 3rd best results are colored with red, green and blue, respectively.
16a
LACASRind
LACASRIE MaanCrpett
LADSONICE
Osjectioon
Orjective
Orsecctin
10
10
10
10
(eccodds)
Time (econds))
(Scconds)
Time (econds)
(a) 'wia'
(b) TDT2'
(c) '20News'
(d) 'sector'
10P
Oloecio
10
EOWC- Brds
Obpectow
10"
Opyeche
5.
10"
10
Time (secorde)
Time (secorda)
Time (secorde)
Time (swconda)
'E2006''
(f) 'MNIST"'
'Gisette"
(h) 'CunCaltech"
Olleciive
Otpect 1e
10
Time Iscconds)
Time (poconds)
'Cidua'
'rando'
Figure 7: The convergence curve for solving L1 regularized SPCA with À = 10.
43

---

Published as a conference paper at ICLR 2026
Obiecive
IIMEC Stst.
10"
10'
Time (secoris)
Time (soconds)
Time (setonds)
Time (octonds)
'wiaa"
(b) 'DT2''
(c) '28News'
'sactor'
10
LABASRind
LADAORICE Miterh
Olgectine
Ostecirs
Oseectre
Ospetirs
10
10"
10
10"
Time (Scoodds)
Time (eecodds)
Time (seconds)
Time Gecoodss
(c) '22006'
(f) 'MNIST
'Gisettee
(h) 'CnnCullechh'
Okyectes
Olpecton
DOEC
Time (secorda)
Time (seconde)
'Cidaa'
'randa'
Figure 8: The convergence curve for solving L, regularized SPCA with A = 50.
Objectiv
Obpective
Oljectos
Obpecivv
10"
12a
1t?
Time (secorde)
Time (secorde)
Time (seconde)
Time (owconde)
'wil''
(b) TDI2'
'20News'
(d) 'sector'
107
10'
Chacch and TEWNE Brds
Olumsion
Obyertow
10"
Time
10"
10"
Time (seconds)
Time (scconds)
Time (sctonds)
(c) 'F2006'
'MNIST"
'Giseme
(h) 'CinCalacch'
Osjectivn
Oxjestion
10
Time (seconds)
Tima (secands)
(i) 'Cifar'
i 'randn'
Figure 9: The convergence curve for solving Legrglllrrzed SPCA with À = 100.
44

---

Published as a conference paper at ICLR 2026
10
Cltmectin
10"
Time (secords)
Time (soconds)
Time (setonds)
Time (occonds)
'wiaa"
(b) 'DT2''
(c) "28News'
'soctor'
10
10"
LABASRind
LADAORICE
Osgectiee
Osgectre
Osjectin
10'
Osgecton
12" Time (scconds)
10"
Time (ecoodas)
Time (ecords))
Time Gecondds
(c) '20006'
(f) 'MNIST
'Gisettee
(h) 'CnnCullechh'
10
Okyectes
Obpecovn
10
Time (seconda)
Time (secorde)
'Cidae'
'randa'
Figure 10: The convergence curve for solving L, regularized SPCA with À = 500.
AUSIM
RADMN
RADMI
10. lime Trrhwo))
tine Srikwl))
W11- 303
02880-1.00 3m-05
T1e-23 20, 0e-00
11e-00
0e-00
15
2077 305
6e-00
3 3100-00 (e-00
01
D0N 5000 1000
2.82--2.30 0--00
14e-30 0m-00
1-33-000
De-m0
3.20--00 01,
26
M0N 5100 1000
-1.3 1.31=43. 0m-00
-3 03, 6e-00
3 3.7ee40 De-th
3. 3.76-00, de
sectredd 4412
4330--33 9m-49
0m-00
30-00
Re-t0
1.2tmm
15-15
64 13: 1000
20 6--00
17--6 (e-m0
1167600 do
Mhiis 60000 TB4
3.57a-02.
0--00
10e-30
0m-00
1.00-00
De-r0
2 777 01,
24
WNIST-600CCS 784
41 130--3 9m-00
48 2b 02, 6m-00
3 347e-0 De-th
3.55-00 2-14
Glerte 1000
41a-1. 0m
0--00
1.80-10
Re-t0
01
Gheze 3000 1000
03 0e-00
6--00
4
(e-m0
3.80-01 76-10
Clae 1000 1000
5p-00
2.91e-30
8e-00
1.50-12
De-00
4 1. 00400
C.n 1080 1000
2 233 9mm00
575
21, 6e-00
12-00, De-th
4. 1.55-0l 06400
mando
-.90-43
2.100-27:
4.76400,
3.24--02
24144
1.000-22
Se-do
Co-10)
4.56-40 (e-40)
2.11e-10
AUSBM
PSS
RADMIX
USCD-R
AUNM
PSM
RADMIX
mmh
(rad)
Trad)
rwd)
[M]
WIG-3199 RC5
Em-10)
Te-4i1
S0,
Io
TD72-500 (000)
3.50p43 00-00
22 8ci00
3.50-40
Be 100
4.00--01
TD72 1000
437
8ci00
455c40
(c-00
4.0co01
2-14
Aacee-113 1000
2.
Ce-te)
Ne-411
30-15
ascred112 1000
Sew50
-1339-50
00 Ce-t0
-1 1396-10
Ne-m
3.Wa-lt S-
20000-2000-1000
4.36d3.000 0o-40
-3.4c-20.00 80i00
Be 100
4.75g-00,
2-14
E2006 2200-2000--0000 1000
3.88243 9c:05
331e-20.8ci000
-13000),
(6-00
4.3ko0) 2-14
MNIST-60000-TA
MNIST-60000-14
S886-02
Glurre sc3c
14-43 56-55
23.58-30,
Co-te)
A--31
5530
Sm-50
-38le-20 Ce-t0
3336-th
Re-m
3.56-tt (b43)
CasCallat- 3000- 1000
4. 00-40
-4.5c-22 8c100
4.74-400
De 100
-.402001
02100
Caycated--2000 1000
4.7,00
9c:05
1L002-20.8cii0
3 (6-00
-Adlcol 0g400
800--300--200
20:100
cundi
56-00
20. Cm-10)
33086-11 S-16
condo
-01a-17 Se-00
37,
1995-10 (e-m
Table 3: Comparisons of objective values and the violation of the nonnegative constraints
(Il mi, X)IlF) for nonnegative PCA for all the compared methods. The 1st, 1 and 3rd best results are colored with red, green and blue, respectively.
45

---

Published as a conference paper at ICLR 2026
Obgedivv
Oltnetive
Olopection
Olteetion
10"
10
10' Time (setondst
Time (seconds)
Time (seconds)
Time (soconds)
(a) "wiaa'
(b) 'DD22'
(c) '20News'
(d) 'sector'
10
Ocyecion
107
Olseccoo
Okty Oceechon
Obgeciis
10"
10-
12
10
10"
10
Time (secorda)
(secorida)
(e) 'E2006''
(f) 'MNIST"
'Gisette"
(h) 'CnnCattcch'
pulppaaad add
Cijectivv
Otpective
10"
Time (pecords)
Time (soconds)
'Cliaa'
"randa"
Figure 11: The convergence curve of the surrogate objective f(X) + 1000|| (0, with X E St(n, r) for solving the nonnegative PCA problem with 10.
Aisoaaad ropc-bidi
Autoihaa copc-boit
Olpecivv
Ollective
OitecTiv
10
10
Time (peconds)
Time (soconds)
Time (soconds)
Time (Sccondst
'Wiaa'
(b) 'TDT2'
'20News'
(d) 'sactor'
10
10
10
Ocjestes
Orgectee
Oojecten
Orgectes
00
10
TOP
10"
10
10'
10a
10
10
10
Time (secords)
Time
Time
(c) 'C2006''
(f) 'MNIST"
'Gisette"
(h) 'CunCallechi
10²
Obyectev
Obpectvv
10*
Time (secords)
Time (secarde)
'Cida''
ii)
'randn'
Figure 12: The convergence curve of the surrogate objective S(X) + 1000 min(0, X)]lF with X e St(n,r) r) for solving the nonnegative PCA problem with 20.
46

---

Published as a conference paper at ICLR 2026
Gooobeed Lop---t d)
10
Oboeddive
Ohtyective
Otopetivon
Olaeetion 12
10+
Time (seconds)
Time (Seconds)
Time (Scconds)
Time isecondst
(a) "wila"
(b) 'DD22'
(c) '20News'
'sector'
107
Ocyecton
Ocyecios
Objectos
6-
10"
10'
10'
l0f
Time (secarda)
Time (secorid))
(e) 'E2006''
(f) 'MNIST
'Gisette"
(h) 'CunCaltcch'
11
Otjectvv
10*
100
Time (seconds)
Time (seconds)
'rando"
Figure 13: The convergence curve of the surrogate objective f(X) + 1000 (0, with X E St(n,r) r) for solving the nonnegative PCA problem with 40.
11
COEC-boit
Olpecive
Ollective
Objecciow 10*
Oltective
10P
1if
Time (seconds)
Time (scconds)
Time (soconds)
Time (sccondst
'wia''
(b) 'TDT2'
'28News'
'sactor'
10
10
Ocjestes
Orgecteen
Ocjecten
Orgectes
00
10'
10
122
10
12²
10' Time (seconds)
(seconds)
Time (secords)
(secords)
(c) 'C2006''
(f) 'MNIST"
'Gisette"
(h) 'CunCallechi
Obyectee
Obpectvve
10*
109
Time (secorde)
Time (secards)
'randn'
Figure 14: The convergence curve of the surrogate objective S(X) + 1000 min(0, X)IlF with X e St(n,r) r) for solving the nonnegative PCA problem with 80.
47


Reviews:
[{"review_id": "ANeic7TAkh", "reviewer": "Reviewer_ufBN", "summary": "The paper proposes a block coordinate descent (OBCD) algorithm for nonsmooth composite optimization on the Stiefel manifold. The method updates $k$ blocks at each iteration while preserving orthogonality, and achieves convergence to a stationary point under suitable assumptions. The authors also design an exact solver for the $k=2$ subproblem when the nonsmooth term $h$ is coordinate-wise separable, by reducing it to a one-dimensional breakpoint-search problem. Theoretical results include global convergence and iteration complexity bounds, and experiments on sparse PCA show promising performance.", "questions": "Please see above.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "* They propose a first BCD method for solving nonsmooth composite optimization problems under orthogonality constraints.\n* The breakpoint-search solver for the $k=2$ subproblem is elegant and provides an exact and efficient solution when $h$ is separable.\n* The experimental results on sparse PCA are convincing.", "weaknesses": "* (Assumption 3 (exact subproblem solution).) This assumption is quite restrictive, as exact solutions are only possible when $h$ is coordinate-wise separable (e.g., $\\ell_0$ or $\\ell_1$ penalties). For general non-separable regularizers such as the nuclear norm, this assumption may not hold. Could the authors discuss (i) whether the framework can be extended to non-separable $h$, and (ii) whether their convergence analysis remains valid if each subproblem is solved only approximately or to a stationary point? A relaxed ``inexact subproblem'' assumption might make the method more generally applicable.\n* (Bounded subgradient assumption (Lemma 4.4).) I think the assumption $||\\partial h(X)||_{sp} \\le \\ell_h$ appears problematic for $h(X)=||X||_0$ (please correct me if I am wrong). However, you use the $\\ell_0$-norm in your experiments. It would be better justify this assumption. \n* (Relation to prior work ([1]).) I think problem (1) is a special case of [1] when $H=L_f \\mathcal{I}_{nr}$. It would strengthen the contribution to explain clearly how this work differs from and improves upon [1].\n* (Experiment) Would the authors also include experimental comparisons with the method proposed in [1]? Such a comparison would better demonstrate the advantages of the proposed OBCD algorithm.\n\n[1] Cheung, Andy Yat-Ming, et al. \"Randomized Submanifold Subgradient Method for Optimization over Stiefel Manifolds.\" arXiv preprint arXiv:2409.01770 (2024).", "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": "XbzKnGmoSL", "reviewer": "Reviewer_TMgr", "summary": "The authors study the problem of minimizing a nonconvex nonsmooth function $F(X)$ over the space of $n \\times r$ orthogonal matrices (the Stiefel manifold $St(n,r)$). This is a challenging class of problems with wide applications. Common approaches for such problems include projection-based methods, Riemannian methods using tangent-space surrogates and retractions, or Block Majorization-Minimization (BMM) methods that iteratively minimize a surrogate directly on the manifold.\n\nThe authors propose an approach, **OBCD**, which falls under the BCD/BMM umbrella but with a novel and distinct design. Instead of updating a manifold block $X_k$ directly, the method subsamples $k$ rows and finds a small $k \\times k$ orthogonal matrix $V$ that *transforms* this block. This \"row-wise\" update ($X^{t+1}(\\mathcal{B},:) \\leftarrow \\overline{V}^{t}X^{t}(\\mathcal{B},:)$) is a key contribution, as it is inherently feasible and avoids the standard tangent-space/retraction machinery.\n\nWhile the general BCD/BMM idea is known, this paper's novelty lies in the specifics of its framework for the Stiefel manifold:\n\n1.  It defines a new optimality condition, the **\"block-k stationary point\" ($BS_k$-point)**.  Theorem 3.6 shows that this condition is **stronger** than the standard critical point condition. The authors justify this by showing their $k=2$ solver uses both rotations and reflections, allowing it to escape suboptimal points.\n\n2.  It provides a **constructive and exact solver** for its nonsmooth subproblem. Appendix B introduces a novel \"Breakpoint Searching Method (BSM)\" that finds the *exact global solution* for the $k=2$ subproblem with $l_0$, $l_1$, or non-negativity regularizers. This is a non-trivial technical result that provides a solid foundation for the algorithm.\n\nOn this theoretical foundation, the authors derive a comprehensive convergence analysis, including an $\\mathcal{O}(1/\\epsilon)$ iteration complexity for an $\\epsilon$-$BS_k$-point (Theorem 4.2) and a full non-ergodic (last-iterate) convergence analysis under the KL property (Theorem 4.10). The presentation is clear, with assumptions formally stated, the algorithm well-defined, and the theoretical claims rigorously established.", "questions": "1. The authors mention in L120 that [Gao et al. 2019] studies a similar problem with columewise updates, whereas the proposed method is rowwise. Besides these algorithmic design choices, what are the important differences? Pros/cons on the computational cost? Theoretical property? \n\n2.  The definition of the \"block-k stationary point\" ($BS_k$-point) (Definition 3.5) is based on $I_k$ being the *global minimizer* of the subproblem. If one only finds a local solution $\\overline{V}^t \\neq I_k$ that satisfies the descent condition (as suggested in Algorithm 1) and the algorithm converges, does the limit point have any meaningful theoretical properties? Does the hierarchy in Theorem 3.6 still hold, or is the $BS_k$-point definition fundamentally tied to the (impractical) global solution of the subproblem for $k > 2$?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "The authors study the problem of minimizing a nonconvex nonsmooth function $F(X)$ over the space of $n \\times r$ orthogonal matrices (the Stiefel manifold $St(n,r)$). This is a challenging class of problems with wide applications. Common approaches for such problems include projection-based methods, Riemannian methods using tangent-space surrogates and retractions, or Block Majorization-Minimization (BMM) methods that iteratively minimize a surrogate directly on the manifold.\n\nThe authors propose an approach, **OBCD**, which falls under the BCD/BMM umbrella but with a novel and distinct design. Instead of updating a manifold block $X_k$ directly, the method subsamples $k$ rows and finds a small $k \\times k$ orthogonal matrix $V$ that *transforms* this block. This \"row-wise\" update ($X^{t+1}(\\mathcal{B},:) \\leftarrow \\overline{V}^{t}X^{t}(\\mathcal{B},:)$) is a key contribution, as it is inherently feasible and avoids the standard tangent-space/retraction machinery.\n\nWhile the general BCD/BMM idea is known, this paper's novelty lies in the specifics of its framework for the Stiefel manifold:\n\n1.  It defines a new optimality condition, the **\"block-k stationary point\" ($BS_k$-point)**. This is a significant contribution, as Theorem 3.6 proves this condition is **provably stronger** than the standard critical point condition. The authors justify this by showing their $k=2$ solver uses both rotations and reflections, allowing it to escape suboptimal points.\n\n2.  It provides a **constructive and exact solver** for its nonsmooth subproblem. Appendix B introduces a novel \"Breakpoint Searching Method (BSM)\" that finds the *exact global solution* for the $k=2$ subproblem with $l_0$, $l_1$, or non-negativity regularizers. This is a non-trivial technical result that provides a solid foundation for the algorithm.\n\nOn this theoretical foundation, the authors derive a comprehensive convergence analysis, including an $\\mathcal{O}(1/\\epsilon)$ iteration complexity for an $\\epsilon$-$BS_k$-point (Theorem 4.2) and a full non-ergodic (last-iterate) convergence analysis under the KL property (Theorem 4.10). The presentation is clear, with assumptions formally stated, the algorithm well-defined, and the theoretical claims rigorously established.", "weaknesses": "### 1. Unclear Practicality for the `k > 2` Case\n\nA primary weakness is the gap between the well-analyzed $k=2$ case and the general $k > 2$ case. The paper's core assumption (Asm-iii) is that the $k \\times k$ nonsmooth subproblem can be solved \"exactly and efficiently.\" The authors provide an impressive, detailed proof of this for $k=2$ using their novel Breakpoint Searching Method (Appendix B).\n\nHowever, for $k > 2$, this assumption is highly questionable. Solving a general $k \\times k$ nonsmooth composite problem over the Stiefel manifold $St(k,k)$ is not a trivial task, and the paper provides no algorithm or justification for it.\n\nThe paper *does* offer a fallback in Algorithm 1, suggesting one can \"Alternatively, find a local solution $\\overline{V}^t$ such that $\\mathcal{K}(\\overline{V}^{t};X^{t},B)\\le\\mathcal{K}(I_{k};X^{t},B)$\". But the paper provides no discussion on *how* to find such a local solution. For a general nonsmooth, nonconvex subproblem, even finding a point that guarantees this simple descent from the identity matrix ($I_k$) is a non-trivial problem in itself. Without a proposed method, the practical application of OBCD for block sizes $k > 2$ remains unclear.\n\n### 2. Overstated Novelty Claim\n\nIn the \"Summary\" of the related work (Section 1.2), the paper claims: \"To our knowledge, this represents the first application of BCD methods to solve nonsmooth composite optimization problems under orthogonality constraints...\". This claim appears to be incorrect. The paper's own literature review (Section 1.2, under \"Minimizing Nonsmooth Functions...\") cites existing work on \"Block Majorization Minimization (BMM) on Riemannian manifolds,\" which directly addresses nonsmooth problems on manifolds, including the Stiefel manifold. This contradiction in the paper's own text weakens the positioning of its contribution. Rather than claiming to be the first of its kind, the authors could emphasize the tailored approach for solving problems on Stiefel manifold efficiently. \n\n### 3. Mismatch Between Theory and Experiments\n\nThere is a disconnect between the theoretical contributions and the experimental validation. The paper provides a strong theoretical justification (in Appendix B) for handling $l_0$, $l_1$, and non-negativity constraints. However, the experiments in Section 5 are conducted *only* on $L_0$-norm-based SPCA. This is a missed opportunity to not validate the new, non-trivial solvers for the $l_1$ and non-negativity constraint problems, which would have made the experimental section much more comprehensive and demonstrated the full power of the proposed subproblem solvers.\n\n### Minor comments\n\nWhile the paper's literature review is adequate, the authors may consider including the following recent paper on Euclidean BMM as part of the BCD section: \n\nHanbaek Lyu and Yuchen Li, “Block majorization-minimization with diminishing radius for constrained nonconvex optimization.”  SIAM Journal on Optimization, Vol. 35, Iss. 2 (2025)", "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": "PrGOAXJwdZ", "reviewer": "Reviewer_Jeff", "summary": "The authors propose an algorithm for minimizing the sum of a smooth and nonsmooth function on the Stiefel manifold.\nRelevant applications include sparse Principal Component Analysis (PCA) with l1 norm and l0 norm regularization.\nThe algorithm is a descent method, and features a block coordinate update, which requires computing the corresponding coordinate block of gradient, as opposed to the full gradient.\nIn theory, the algorithm converges in expectation to a relevant point with classical complexity (Th. 4.2), and under additional assumptions (including continuity of $h$, which rules the important l0 norm example and leaves l1 norm) the last iterate converges to a relevant point with classical complexity. We note that the notion of \"relevant point\" is slightly stronger than the usual notion of critical point (0 in limiting subdifferential).\nIn practice, the algorithm performs well on PCA with l0 regularization.", "questions": "Any comment and detail on each of the listed weak points is welcome.\n\nBelow are some suggestions in the form of questions; any answer on the following questions is welcome, but no answer is also fine.\n- connection to (variable metric) proximal gradient method? That would connect assumption (iii) to the notion of \"prox-friendly\" nonsmooth function $h$, standard in nonsmooth optimization. That may also help and connect\n- Th. 3.6: globally optimal points are BS2-points. Are they also BSk points, for $2 \\le k \\le n-1$?\n- Assumption (iii) is reminiscent of the  might appear more natural if connected to the proximal gradient operator.\n- optimization methods on problems with l0-norm regularization usually face the combinatorial difficulty of the l0 norm in some way. I am surprised that it doesn't show in your analysis. Do you have some intuition to share on this? I am also surprised that the complexity of OBCD with parameter $k$ in Theorem 4.2 does not depend on $k$. Again, do you have some intuition to share on this?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "- the submission is clearly written,\n- the problem is motivated by a relevant application (sparse PCA with l0 and l1 norm regularization),\n- the algorithm is derived and formulated in a rigorous and clear way (up to one detail, discussed below),\n- the theoretical analysis presents a notion of stationary points (BS$k$-points) that is stronger than the customary critical point (zero in subdifferential), and a proof that the algorithm converges to such points on average. I am quite familiar with composite optimization, but less so on block coordinate schemes so this idea may be common in that field.\n- the algorithm combines known ideas (block coordinate descent, majorization minimization) suitably, and successfully applies these ideas to minimization with l0 norm regularization. As far as I know, this task is highly challenging, and few methods for it are available even in the unconstrained case.", "weaknesses": "Here is a list of points that prevent me from providing a better assessment to the submission.\n- The nonsmooth function $h$ has strong, restrictive assumptions. This severely limits the applicability of the method.\n  + Assumption (ii): $h$ is assumed to be coordinate separable, with the same expression for each coordinate.\n  + Assumption (iii) on function $h$ is shown to be tractable in theory for three values of $h$ only (the $\\ell1$ norm, the indicator of a polyhedron, and the $\\ell0$ norm when $k=2$).\n- l. 122-128: I disagree with part of the positioning relative with the literature:\n  + contrary to what point (ii) implies, I think that the proposed method is not applicable to general nonsmooth composite problems. Indeed it relies on an assumption on the nonsmooth function, which is arguably restrictive (see point below) but certainly not valid for \"general nonsmooth composite problems\".\n  + point (iii) asserts that it is a \"limitation\" for methods to be \"infeasible\", that is to generate iterates that are only asymptotically feasible. Yet, infeasible methods are accepted for solving large-scale optimization problems under general constraints (e.g. with the Augmented Lagrangian method), and Stiefel manifold constraints specifically (e.g. with the recent Landing methods, see next point).\n  + point (iv), \"they often lack rigorous convergence guarantees\" should be more precise on the type of convergence guarantees that previous methods fail (for instance: convergence of average, last iterate, complexity guarantee).\n- l. 78-107: Literature review misses one recent line of work of so-called \"landing methods\" for optimization on Stiefel manifolds: [1-4]. These methods also target large-scale minimization on the Stiefel manifold; the paper should discuss these methods.\n- Algorithm 1: The algorithm seems to depend on a parameter $\\alpha$ (l. 194, Lemma 2.3, Th. 4.2, among other occurrences), yet this parameter does not appear in the algorithm statement, and it is unclear what condition $\\alpha$ should meet so that the convergence guarantees hold, and particularly so the sufficient decrease condition of Theorem 4.2(a).\n- l. 374: the condition $\\| \\partial h(X) \\|_{sp}<l_h$ is ill-posed: the \"sp\" norm of a set is not defined. Besides, this condition may appear more natural if it were implied by Lipschitz continuity of $h$.\n- l. 392: the assumption that $F_i$ is a KL function is not discussed at all. It is thus unclear whether it holds for any of the application cases, and the three discussed $h$ functions.\n- l. 395: Proposition 4.8, from previous work, is stated without providing a precise reference, that is, both the paper and statement references.\n- l. 404 & 412: \"the continuity assumption made in lemma 4.4\" is unclear. To what condition does this sentence refer to exactly? Maybe writing an additional assumption for this would clarify the situation.\n- Experiments: it is unclear whether the three baselines generate feasible iterates. In addition, these three baselines involve operator splitting. Such methods usually involve a reformulation of the problem with additional constraints, which are satisfied asymptotically only. Finally, it is not clear whether the three baselines have convergence guarantees on  the problem with $h = \\| \\cdot \\|_{0}$, which is discontinuous; this aspect conditions the interpretation of the whole experimental section.\n\nReferences:\n- [1] Goyens, Absil & Feppon (2026) Geometric Design of the Tangent Term in Landing Algorithms for Orthogonality Constraints, Springer Nature Switzerland.\n- [2] Ablin & Peyré (2022) Fast and Accurate Optimization on the Orthogonal Manifold without Retraction, PMLR.\n- [3] Vary, Ablin & Gao et al. (2024) Optimization without Retraction on the Random Generalized Stiefel Manifold, PMLR.\n- [4] Gao, Vary & Ablin et al. (2022) Optimization Flows Landing on the Stiefel Manifold$\\star$, IFAC-PapersOnLine.\n\n\nMinor points, that do not impact the assessment of the paper:\n- l. 35: clarity would improve by mentioning explicitly that $f$ is assumed to be differentiable and $H$-smooth\n- l. 39, footnote: Is $f$ assumed to have expression $1/2 \\|X\\|_{H}^{2}$ on the whole paper?\n- l. 40 & 296: what are the definitions of $F$ \"closed\", and $F$ lower semicontinuous? How do they differ?\n- l. 230: is the\n- l. 305: The writing of Definition 3.3 is confusing: the sentence \"Furthermore $\\lambda \\in [\\partial F(X)]^\\top X$\" can read as an additional condition for a point $X$ to be critical, but the notion  of criticality does not involve $\\lambda$.\n- l. 313-318: what is the purpose of this paragraph?", "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": "Me83E95duK", "reviewer": "Reviewer_X2XH", "summary": "A new approach called OBCD is proposed for minimizing the sum of a smooth and a nonsmooth convex functions over the nonconvex set of orthogonal matrices.", "questions": "N/A", "limitations": "", "rating": 8, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Optimization under orthogonality constraints is an important topic with many applications. The paper is well written. My expertise of the topic is limited, though. I know about projection methods, as there is a surge of interest for orthogonalization recently in the context of the Muon optimizer, see for instance Grishina et al. \"Accelerating Newton-Schulz Iteration for Orthogonalization via Chebyshev-type Polynomials\". The proposed method maintains feasibility throughout the process, which is a different approach than relying on approximate orthogonalization.", "weaknesses": "I don't see limitations.", "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)