Block Coordinate Descent Methods for Optimization under J-Orthogonality Constraints with Applications

Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 07:07:58

Abstract: The J-orthogonal matrix, also referred to as the hyperbolic orthogonal matrix, is a class of special orthogonal matrix in hyperbolic space, notable for its advantageous properties. These matrices are integral to optimization under J-orthogonal constraints, which have widespread applications in statistical learning and data science. However, addressing these problems is generally challenging due to their non-convex nature and the computational intensity of the constraints. Currently, algorithms for tackling these challenges are limited. This paper introduces \textbf{JOBCD}, a novel Block Coordinate Descent method designed to address optimizations with J-orthogonality constraints. We explore two specific variants of \textbf{JOBCD}: one based on a Gauss-Seidel strategy (\textbf{GS-JOBCD}), the other on a variance-reduced and Jacobi strategy (\textbf{VR-J-JOBCD}). Notably, leveraging the parallel framework of a Jacobi strategy, \textbf{VR-J-JOBCD} integrates variance reduction techniques to decrease oracle complexity in the minimization of finite-sum functions. For both \textbf{GS-JOBCD} and \textbf{VR-J-JOBCD}, we establish the oracle complexity under mild conditions and strong limit-point convergence results under the Kurdyka-Lojasiewicz inequality. To demonstrate the effectiveness of our method, we conduct experiments on hyperbolic eigenvalue problems, hyperbolic structural probe problems, and the ultrahyperbolic knowledge graph embedding problem. Extensive experiments using both real-world and synthetic data demonstrate that \textbf{JOBCD} consistently outperforms state-of-the-art solutions, by large margins.

Main paper: 
 introduction: 
 A matrix X ∈ R n×n is a J-orthogonal matrix if X T JX = J, where J = [ Ip 0 0 -In-p ], and I p is a p × p identity matrix. Here, J ∈ R n×n is the signature matrix with signature (p, n -p). In this paper, we mainly focus on the following optimization problem under J-orthogonality constraints:EQUATIONHere, f (X) could have a finite-sum structure, each component function f i (X) is assumed to be differentiable, and N is the number of data points. For brevity, the J-orthogonality constraint 1) is rewritten as X ∈ J .X T JX = J in Problem (We impose the following assumptions on Problem (1) throughout this paper. (A-i) For any matrices X and X + , we assume f i : R n×n → R is continuously differentiable for some symmetric positive semidefinite matrix H ∈ R nn×nn that:EQUATIONfor all i ∈ [N ], where ∥H∥ ≤ L f for some constant L f > 0 and ∥X∥ 2 H ≜ vec(X) ⊤ H vec(X).This further implies that: ∥∇f i (X) -∇f i (X + )∥ F ≤ L f ∥X -X + ∥ F for all i ∈ [N ]. Importantly, the function f (X) = 1 2 tr(X ⊤ CXD) = ∀X, X + , f (X + ) = Q(X + ; X) in (2) , where C ∈ R n×n and D ∈ R n×n are arbitrary symmetric matrices. (A-ii) The function f i (X) is coercive for all i ∈ N , that is, lim ∥X∥ F →∞ f i (X) = ∞, ∀i.Problem (1) defines an optimization framework that is fundamental to a wide range of models in statistical learning and data science, including hyperbolic eigenvalue problem [6, 43, 40] , hyperbolic structural probe problem [20, 7] , and ultrahyperbolic knowledge graph embedding [48] . Additionally, it is closely related to machine learning in hyperbolic spaces, including Lorentz model learning [35, 50, 8] and ultrahyperbolic neural networks [27, 54, 42] . It also intersects with hyperbolic linear algebra [3, 21] , addressing problems such as the indefinite least squares problem, hyperbolic QR factorization, and indefinite polar decomposition.▶ Block Coordinate Descent Methods. Block Coordinate Descent (BCD) is a well-established iterative algorithm that sequentially minimizes along block coordinate directions. Its simplicity and efficiency have led to its widespread adoption in structured convex applications [37] . Recently, BCD has gained traction in non-convex problems due to its robust optimality guarantees and/or excellent empirical performance in areas including optimal transport [22] , matrix optimization [12] ,fractional minimization [52] , deep neural networks [5, 53, 32] , federated learning [47] , black-box optimization [4] , and optimization with orthogonality constraints [51, 14] . To our knowledge, this is the first application of BCD methods to optimization under J-orthoginality constraints, with a focus on analyzing their theoretical guarantees and empirical efficacy.▶ Minimizing Smooth Functions under J-Orthogonality Constraints. The J-orthogonal matrix belongs to a subset of generalized orthogonal matrices [16, 36, 23] . However, projecting onto the J-orthogonality constraint poses challenges, complicating the extension of conventional optimization algorithms to address optimization problems under these constraints [1, 16] . This contrasts with computing orthogonal projections using methods such as polar or SVD decomposition, or approximating them via QR factorization. Existing methods for addressing Problem (1) can be categorized into three classes. (i) CS-Decomposition Based Methods. These approaches involve parameterizing four orthogonal matrices (as described in Proposition 2.2) and subsequently minimizing a smooth function over these matrices in an alternating fashion. The involvement of 3 × 3 block matrices makes the implementation of these methods very challenging. Consequently, the work of [48] focuses on optimizing a reduced subspace of the CS decomposition parameters, albeit at the expense of losing some degrees of freedom. (ii) Unconstrained Multiplier Correction Methods [31, 13, 14] . These methods leverage the symmetry and explicit closed-form expression of the Lagrangian multiplier at the first-order optimality condition. Consequently, they address an unconstrained problem, resulting in efficient first-order infeasible approaches. (iii) Alternating Direction Method of Multipliers [19] .This method reformulates the original problem into a bilinear constrained optimization problem by introducing auxiliary variables. It employs dual variables to handle bilinear constraints, iteratively optimizing primal variables while keeping other primal and dual variables fixed, and using a gradient ascent strategy to update the dual variables. This approach has become widely adopted for solving general nonconvex and nonsmooth composite optimization problems. Notably, all the aforementioned methods solely identify critical points of Problem (1) .▶ Finite-Sum Problems via Stochastic Gradient Descent. The finite-sum structure is prevalent in machine learning and statistical modeling, facilitating decomposition into smaller, more manageable components. This property is advantageous for developing efficient algorithms for large-scale problems, such as Stochastic Gradient Descent (SGD). Reducing variance is crucial in SGD because it can lead to more stable and faster convergence. Various techniques, such as mini-batch SGD, momentum methods, and variance reduction methods like SAGA [10] , SVRG [25] , SARAH [34] , SPIDER [11, 44] , SNVRG [55] , and PAGE [30] , have been developed to address this issue. Additionally, SGD for minimizing composite functions has also been investigated by the authors [15, 24, 29] .This paper makes the following contributions. (i) Algorithmically: We introduce the JOBCD algorithm, a novel Block Coordinate Descent method specifically designed to tackle optimizations constrained by J-orthogonality. We explore two specific variants of JOBCD, one based on a Gauss-Seidel strategy (GS-JOBCD), the other on a variance-reduced and Jacobi strategy (VR-J-JOBCD). Notably, VR-J-JOBCD incorporates a variance-reduction technique into a parallel framework to reduce oracle complexity in the minimization of finite-sum functions (See Section 2). (ii) Theoretically: We provide comprehensive optimality and convergence analyses for both Algorithm 1: GS-JOBCD: Block Coordinate Descent Methods using a Gauss-Seidel Strategy for Solving Problem (1) Init.: Set X 0 to satisfy J-orthogonality constraints (e.g., via Hyperbolic CS Decomposition). for t from 0 to T do (S1) Choose a coordinate B t with |B t | = 2 from the set Ω randomly and uniformly (with replacement) for the t-th iteration. Denote B = B t . (S2) Choose a matrix Q ∈ R 4×4 using Formula (4). (S3) Solve the following small-size subproblem globally.EQUATION= arg minEQUATIONwhereP ≜ [∇f (X t )(X t ) T ] BB -mat( Q vec(I 2 )), Q = Q + θI and c ≜ f (X t ) -⟨I 2 , [∇f (X t )(X t ) T ] BB ⟩ + 1 2 ∥I∥ 2 Q is a constant. (S4) X t+1 (B, :) = V t X t (B, :)end Proposition 2.2. (Hyperbolic CS Decomposition [41] ) Let V be J-orthogonal with signature (p, n -p). Assume that n -p ≤ p. Then there exist vectors ċ, ṡ ∈ R n-p with ċ ⊙ ċṡ ⊙ ṡ = 1, and orthogonal matrices U 1 , V 1 ∈ R p×p and U 2 , V 2 ∈ R (n-p)×(n-p) such that:V = [ U 1 0 0 U 2 ][ Diag( ċ) 0 Diag( ṡ) 0 I p-(n-p) 0 Diag( ṡ) 0 Diag( ċ) ][ V T 1 0 0 V T 2 ].Applying Proposition 2.2 with n = 2, p = 1, andU 1 = U 2 = V 1 = V 2 = ±1, c2 -s2 = 1 with c, s ∈ R, we parametrize V as: V = ( ±1 0 0 ±1 ) • ( c s s c ) • ( ±1 0 0 ±1), where we denote s as sinh(µ), c as cosh(µ), and t as tanh(µ)for some µ ∈ R, for simplicity of notation. It is not difficult to show that Problem (6) reduces to the following one-dimensional search problem: μ ∈ minWe apply a breakpoint search method to solve Problem (7) . For simplicity, we provide an analysis only for the first case. A detailed discussion of all four cases can be found in Appendix Section B.1.For the case where V = ( c s s c ), Problem (7) reduces to the following problem: min c,s a c +b s +c c2 +d c s +e s2 ,where a = P 11 + P 22 , b = P 12 + P 21 , c = 1 2 ( Q11 + Q41 + Q14 + Q44 ), d = 1 2 ( Q21 + Q31 + Q12 + Q42 + Q13 + Q43 + Q24 + Q34 ), and e = 1 2 ( Q22 + Q32 + Q23 + Q33 ). Then we perform a substitution to convert Problem (8) into an equivalent problem that depends on the trigonometric functions:(i) c2 = 1 1- t2 ; (ii) s2 = t2 1- t2 ; (iii) t = s c .The following lemma provides a characterization of the global optimal solution for Problem (8) .Lemma 2.3. (Proof in Section C.2) We let F (c, s) ≜ ac + bs + cc 2 + dcs + es 2 . The optimal solution μ to Problem (8) can be computed as:[cosh(μ), sinh(μ)] ∈ arg min [c,s] F (c, s), s. t. [c, s] ∈ {[ 1 √ 1-( t+) 2 , t+ √ 1-( t+) 2 ], [ -1 √ 1-( t-) 2 , -t- √ 1-( t-) 2 ]}, where t+ ∈ arg min t p(t) ≜ a+bt √ 1-t 2 + w+dt 1-t 2 ; t-∈ arg min t p(t) ≜ -a-bt √ 1-t 2 + w+dt 1-t 2 .Here w = c + e.We now describe how to find the optimal solution t+ , where t+ ∈ arg min t p(t) ≜ a+bt √ 1-t 2 + w+dt 1-t 2 ; this strategy can naturally be extended to find t-. Initially, we have the following first-order optimality conditions for the problem:0 = ∇p(t) =[b(1 -t 2 ) + (a + bt)t] √ 1 -t 2 + [d(1 -t 2 ) + (w + dt)(2t)] ⇔ dt 2 + 2wt + d = -[b + at] √ 1 -t 2 .Squaring both sides yields the following quartic equation:c 4 t 4 + c 3 t 3 + c 2 t 2 + c 1 t + c 0 = 0, where c 4 = d 2 + a 2 , c 3 = 4wd + 2ab, c 2 = 4w 2 + 2d 2 -a 2 + b 2 , c 1 = 4wd -2ab, c 0 = d 2 -b 2. This equation can be solved analytically by Lodovico Ferrari's method [46] , resulting in all its real roots { t1 , t2 , . . . , tj } with 1 ≤ j ≤ 4.For the second and third cases, Problem (6) essentially boils down to optimization under orthogonality constraints. The work of [51] derives a breakpoint search method for finding the optimal solution for Problem (6) with J BB ∈ {(1 0 0 1 ), ( -1 0 0 -1 )} using the Givens rotation and Jacobi reflection matrices. 
variance-reduced jacobi block coordinate descent algorithm: 
 This subsection proposes the VR-J-JOBCD algorithm, a randomized block coordinate descent method derived from GS-JOBCD. Importantly, by leveraging the parallel framework of a Jacobi strategy [17, 9] , VR-J-JOBCD integrates variance reduction techniques [39, 30, 18] to decrease oracle complexity in the minimization of finite-sum functions. This makes the algorithm effective for minimizing large-scale problems under J-orthogonality constraints.Notations. We assume n is an even number in this paper. We create (n/2) pairs by non-overlapping grouping of the numbers in any arbitrary combination, with each pair containing two distinct numbers from the set [n]. It is not hard to verify that such grouping yieldsC J = (n!)/(2 n/2 n 2 !) possiblecombinations. The set of these combinations is denoted asΥ ≜ { Bi } C J i=1 ≜ { B1 , B2 , . . . , BC J } 1 .▶ Variance Reduction Strategy. We incorporate state-of-the-art variance reduction strategies from the literature [30, 5] into our algorithm to solve Problem (1). These methods iteratively generate a stochastic gradient estimator as follows:Gt = 1 b i∈S t + ∇f i (X t ), with probability p; Gt-1 + 1 b ′ i∈S t 1 b i∈S 0 + ∇f i (X 0). We drop the superscript t for {S t + , S t * } as t can be inferred from context. We only focus on the default setting that [30, 5] : b = N , b ′ = √ b and p = b ′ b+b ′ .▶ Jacobi Block Coordinate Descent Method. The proposed algorithm is built upon the parallel framework of a Jacobi strategy. In each iteration t, we randomly and uniformly (with replacement)select a coordinate set B t ≜ {B t (1) , B t (2) , • • • , B t ( n 2 ) } from the set Υ with B t ∈ N n 2 ×2 and B t (i) ∈ N 2 .For all t, we have:B t (i) ∩ B t (j) = ∅ and ∪ n/2 i=1 (B t (i) ) = [n].We drop the superscript t if t can be inferred from context.The following lemma shows how to choose a suitable matrix Q so that the Jacobi strategy can be applied.Lemma 2.4. (Proof in Section C.3) We let B t ≜ {B t (1) , B t (2) , • • • , B t ( n 2 ) } ∈ Υ for all t. We let Q = ςI 4 ,where ς is some suitable constant with ς ≤ L f . For any B t (i) and B t (j) with i ̸ = j, their corresponding objective functions as in Equation (3) are independent.We consider the following block coordinate update rule in VR-J-JOBCD:X t+1 ⇐ X t B (V : ) ≜ X t + [ n/2 i=1 U B (i) (V i -I 2 )U ⊤ B (i) ]X t .The following lemma provides properties of this rule.Lemma 2.5. (Proof in Section C.4) We let B ∈ Υ, V i ∈ J B (i) , X ∈ J , and i ∈ [ n 2 ]. We define X + ≜ XB (V : ) ≜ X + [ n/2 i=1 U B (i) (V i -I 2 )U ⊤ B (i) ]X. We have: (a) n 2 i=1 ∥U B (i) (V i -I 2 )U ⊤ B (i) X∥ 2 F = ∥ n 2 i=1 U B (i) (V i -I 2 )U ⊤ B (i) X∥ 2 F . (b) ∥X + -X∥ 2 F ≤ ∥X∥ 2 F • n/2 i=1 ∥V i -I 2 ∥ 2 F . (c) ∥X + -X∥ 2 H ≤ n/2 i=1 ∥V i -I 2 ∥ 2 Q with Q = ςI 4 . (d) For all G ∈ R n×n , it follows that: 2 n/2 i=1 ⟨V i -I 2 , [(∇f (X) -G)X ⊤ ] B (i) B (i) ⟩ ≤ ∥X∥ 2 F n/2 i=1 ∥V i -I 2 ∥ 2 F + ∥[∇f (X) - G]∥ 2 F .▶ The Main Algorithm. Using the update rule above, we consider the following iterative procedure:X t+1 ⇐ X t B (V : ), where Vt : ∈ arg min V: f ( X t B (V : )).We establish the majorization function for f ( X t B (V : )), as follows:f ( X t B (V : )) ① ≤ f (X t ) + ⟨ X t B (V : ) -X t , ∇f (X t )⟩ + 1 2 ∥ X t B (V : ) -X t ∥ 2 H ② ≤ f (X t ) + n/2 i=1 {⟨V i -I 2 , [∇f (X)(X) ⊤ ] B (i) B (i) ⟩ + 1 2 ∥V i -I 2 ∥ 2 (θ+ς)I } (10)where step ① uses the results of telescoping Inequality (2) over i from 1 to N ; step ② uses X t+1 -X t = [ n/2 i=1 U B (i) (V i -I 2 )U ⊤ B (i) ]X t , Claim (c) of Lemma 2.5, θ ≥ 0, and Q = ςI.Instead of computing the exact Euclidean gradient ∇f (X t ) as GS-JOBCD, VR-J-JOBCD maintains and updates a recursive gradient estimator Gt using a variance-reduced strategy as in Formula (9).We consider minimizing the following function instead of the one on the right-hand side of Inequality (10) :EQUATION)Here, T (V : ; X t , B t ) can be termed as a stochastic majorization function of f ( X t B (V : )) at the current solution X t . Therefore, we can consider the following optimization problem to find {V : } using:Vt : ∈ arg min V: T (V : ; X t , B t), which can be decomposed into (n/2) independent subproblems and solved in parallel. It is important to note that each V i in Problem ( 12) is identical to Problem (6), which can be efficiently solved in O(1) using the breakpoint search method, as in GS-JOBCD.We summarize the proposed VR-J-JOBCD in Algorithm 2. Notably, when N = 1, VR-J-JOBCD simplifies to a direct Jacobi strategy for solving Problem (1), which we refer to as J-JOBCD.Algorithm 2: VR-J-JOBCD: Block Coordinate Descent Methods using a variance-reduced and Jacobi strategy for Solving Problem 1 Init.: Set X 0 to satisfy J-orthogonality constraints (e.g., via Hyperbolic CS Decomposition). for t from 0 to T do (S1) Choose a coordinate B t from the set Υ randomly and uniformly (with replacement) for the t-th iteration. Denote B = B t . In our implementation, we simply randomly permute the set {1, 2, ..., n} and then output the grouping {[1, 2], [3, 4] , [5, 6] , for i = 1 to n/2 in parallel doEQUATIONEQUATIONwhereP i ≜ [∇f (X t )(X t ) T ] B (i) B (i) -mat( Q vec(I 2 )) -θI 2 , Q = (ζ + θ)I. (S4)Update the solution X t+1 in parallel as follows:for i = 1 to n/2 in parallel do X t+1 (B (i) , :) = Vt i X t (B (i) , :) endQuestion: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes] Justification: Please refer to the assumptions made for the optimization problem outlined in the introduction and Section 4.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.• The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.• The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.• The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.• The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.• If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.• While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
optimality analysis: 
 This section provides an optimality analysis for the proposed algorithms.Initially, we define the first-order optimality condition for Problem (1) . Since the matrix X T JX is symmetric, the Lagrangian multiplier Λ corresponding to the constraints X T JX = J is also asymmetric matrix. The Lagrangian function of problem (1) is L(X, Λ) = f (X) -1 2 ⟨Λ, X T JX -J⟩.We obtain the following lemma for the first-order optimality condition for Problem (1). 0 = ∇ J f ( X) ≜ ∇f ( X) -J X[∇f ( X)] ⊤ XJ.The associated Lagrangian multiplier can be computed asΛ = J XT ∇f ( X). (b)The critical point condition is equivalent to the requirement that the matrix X∇f ( X) T J is symmetric, which is expressed asXG T J = [XG T J] T .Remarks. While our results in Lemma 3.1 show similarities to existing works focusing on problems under orthogonality constraints [45] , this study marks the first investigation into the first-order optimality condition for optimization problems under J-orthogonality constraints.The following definition is useful in our subsequent analysis of the proposed algorithms. I 2 ∈ arg min V∈JB G(V; Ẍ, B).The following theorem shows the relation between critical points and BS-points.Theorem 3.3. (Proof in Section D.2) Any BS-point is a critical point, while the reverse is not necessarily true.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?• The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).• The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)• The assumptions made should be given (e.g., Normally distributed errors).• It should be clear whether the error bar is the standard deviation or the standard error of the mean.• It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified.• For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates).• If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
convergence analysis: 
 This section provides a convergence analysis for GS-JOBCD and VR-J-JOBCD.For GS-JOBCD, the randomness of output (V t , X t+1 ) for all t are influenced by the random variableξ t ≜ (B 1 ; B 2 ; • • • ; B t ).For VR-J-JOBCD, the randomness of output ( Vt : , X t+1 ) are influenced by the random variablesι t ≜ (B 1 , S 1 + , S 1 * ; B 2 , S 2 + , S 2 * ; • • • ; B t , S t + , S t * ).We denote X as the global optimal solution of Problem (1). To simplify notations, we define:u t = ∥ Gt -∇f (X t )∥ 2 F , and ∆ i = f (X i ) -f ( X).We impose the following additional assumptions on the proposed algorithms.Assumption 4.1. There exists constants {X, V} that: ∥X t ∥ F ≤ X, and ∥V t ∥ F ≤ V for all t.Assumption 4.2. There exists a constant G that: ∥∇f (X t )∥ F ≤ G, and ∥ Gt ∥ F ≤ G for all t.Assumption 4.3. For anyX ∈ R n×n , E i [∥∇f i (X t ) -∇f (X t )∥ 2 F ] ≤ σ 2, where i is drawn uniformly at random from [N ].Remarks. (i) Assumption 4.1 is satisfied as the function f i (X) is coercive for all i. (ii) Assumption 4.2 imposes a bound on the (stochastic) gradient, a fairly moderate condition frequently employed in nonconvex optimization [26] . (iii) Assumption 4.3 ensures that the variance of the stochastic gradient is bounded, which is a common requirement in stochastic optimization [30, 5] .We define the ϵ-BS-point as follows.Definition 4.4. (ϵ-BS-point) Given any constant ϵ > 0, a point Ẍ is called an ϵ-BS-point if: E( Ẍ) ≤ ϵ.Here, E(X) is defined asE(X) ≜ 1 C 2 n C 2 n i=1 dist(I 2 , arg min V G(V; X, B i )) 2 for For GS-JOBCD, while it is defined as E(X) ≜ 1 C J C J i=1 E ι t [dist(I 2 , arg min V: T (V : ; X, Bi )) 2 ] for VR-J-JOBCD,where the expectation is with respect to the randomness inherent in the algorithm [30] .We have the following useful lemma for VR-J-JOBCD. E ι t [u k ] of the gradient estimators { Gt } of Algorithm 2 is bounded by: E ι t [u t ] ≤ p(N -b) b(N -1) σ 2 + (1 -p)E ι t-1 [u t-1 ]+ L 2 f X 2 (1-p) b ′ E ι t-1 [ n/2 i=1 ∥V t-1 i -I 2 ∥ 2 F ]The following two theorems establish the iteration complexity (or oracle complexity) for GS-JOBCD and VR-J-JOBCD. Remark. Theorems 4.6 and 4.7 demonstrate that the arithmetic operation complexity of GS-JOBCD is linearly dependent on N , while VR-J-JOBCD is linearly dependent on √ N . Therefore, VR-J-JOBCD reduces the iteration complexity significantly.We prove algorithms achieve strong convergence based on a non-convex analysis tool called Kurdyka-Łojasiewicz inequality [2] .We impose the following assumption on Problem (1) .Assumption 4.8. (Kurdyka-Łojasiewicz Property). Assume that f • (X) = f (X) + I J (X) is a KL function.For all X ∈ dom f • , there exists σ ∈ [0, 1), η ∈ (0, +∞] a neighborhood Υ of X and a concave and continuous function φ(t) = ct 1-σ , c > 0, t ∈ [0, η) such that for all X ′ ∈ Υ and satisfiesf • (X ′ ) ∈ (f • (X), f • (X) + η), the following holds: dist(0, ∇f • (X ′ ))φ ′ (f • (X ′ ) -f • (X)) ≥ 1.We establish strong limit-point convergence for VR-J-JOBCD and GS-JOBCD.Theorem 4.9. (Proof in Section E.5, a Finite Length Property). The sequence {X t } ∞ t=0 of GS-JOBCD has finite length property that: ∀t, t i=1 E ξ t [∥X t+1 -X t ∥ F ] ≤ O(φ(∆ 1 )) < +∞,∀t, t i=1 E ι t [∥X t+1 -X t ∥ F ] ≤ O( φ(∆1) N 1/4) < +∞, where φ(•) is the desingularization function defined in Assumption 4.8. 
applications and numerical experiments: 
 This section demonstrates the effectiveness and efficiency of JOBCD on three optimization tasks:(i) the hyperbolic eigenvalue problem, (ii) structural probe problem, and (iii) Ultra-hyperbolic Knowledge Graph Embedding problem. We provide experiments for the last problem in Section F.2.▶ Application to the Hyperbolic Eigenvalue Problem (HEVP). The hyperbolic eigenvalue problem refers to the generalized eigenvalue problem in hyperbolic spaces [40] . This problem is a fundamental component in machine learning models, such as Hyperbolic PCA [43, 6] . Given a data matrix D ∈ R m×n and a signature matrix J with signature (p, n -p), HEVP can be formulated as the following optimization problem: min X -tr(X T D T DX), s. t. X T JX = J.▶ Application to the Hyperbolic Structural Probe Problem (HSPP). The Structure Probe (SP) is a metric learning model aimed at understanding the intrinsic semantic information of large language models [20] [7] . Given a data matrix D ∈ R m×n and its associated Euclidean distance metric matrix T ∈ R m×m , HSPP employs a smooth homeomorphic mapping function φ(•) to project the data D into ultra-hyperbolic space. Subsequently, it seeks an appropriate linear transformation X ∈ R n×n constrained within a specific structure X ∈ J , such that the resulting transformed data Q ≜ φ(D)X ∈ R m×n exhibits similarity to the original distance metric matrix T under the ultra-hyperbolic geodesic distanced α (Q i: , Q j: ), expressed as T i,j ≈ d α (Q i: , Q j: ) for all i, j ∈ [m],where Q i: is i-th row of the matrix Q ∈ R m×n . This can be formulated as the following optimization ▶ Datasets. To generate the matrix D ∈ R m×n , we use 8 real-world or synthetic data sets for both HEVP and HSPP tasks: 'Cifar', 'CnnCaltech', 'Gisette', 'Mnist', 'randn', 'Sector', 'TDT2', 'w1a'.problem: min X 1 m 2 i,j∈m (T i,j -d α (Q i: , Q j: )) 2 , s. t. Q ≜ φ(D)X, X ∈ J .We randomly extract a subset from the original data sets for the experiments.▶ Compared Methods. We compare GS-JOBCD and VR-J-JOBCD with 3 state-of-the-art optimization algorithms under J-orthogonality constraints. (i) The CS Decomposition Method (CSDM) [48] . (ii) Stardard ADMM (ADMM) [19] . UMCM: Unconstrained Multiplier Correction Method [31, 13] . ▶ Experiment Settings. All methods are implemented using Pytorch on an Intel 2.6 GHz processor with an A40 (48GB). For HSPP, we fix α to 1. Each method employs the same random J-orthogonal matrix. The built-in solver Admm is used to solve the unconstrained minimization problem in CSDM.We provide our code in the supplemental material. ▶ Experiment Results. Table 1 and Figure 1 display the accuracy and computational efficiency for HEVP, while Figure 2 presents the results for HSPP, leading to the following observations: (i)GS-JOBCD and JJOBCD consistently deliver better performance than the other methods. (ii) Other methods frequently encounter poor local minima, whereas GS-JOBCD effectively escapes these minima and typically achieves lower objective values, aligning with our theory that our methods locate stronger stationary points. (iii) VR-J-JOBCD outperforms both J-JOBCD and CSDM when dealing with a large dataset characterized by an infinite-sum structure. 
conclusions: 
 In this paper, we propose a new approach JOBCD, which is based on block coordinate descent, for solving the optimization problem under J-orthogonality constraints. We discuss two specific variants of JOBCD: one based on a Gauss-Seidel strategy (GS-JOBCD), the other on a variance-reduced Jacobi strategy. Both algorithms capitalize on specific structural characteristics of the constraints to converge to more favorable stationary solutions. Notably, VR-J-JOBCD incorporates a variancereduction technique into a parallel framework to reduce oracle complexity in the minimization of finite-sum functions. For both GS-JOBCD and VR-J-JOBCD, we establish the oracle complexity under mild conditions and strong limit-point convergence results under the Kurdyka-Lojasiewicz inequality. Some experiments on the hyperbolic eigenvalue problem and structural probe problem show the efficiency and efficacy of the proposed methods. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes] Justification:Guidelines: We have outlined the types of compute workers, detailing CPU and memory specifications.• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.• The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.• The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines? 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [NA]Justification: The paper addresses theoretical questions on algorithm complexity, which, to the best of our knowledge, pose no negative social impact.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.• Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.• The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.• The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.• If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: The paper poses no such risks.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.• Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.• We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort. Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL.• The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. 
 
 Appendix: Lemma A.1. (Lemma 6.6 of [51] ) For any W ∈ R n×n , we have:C k n i=1 ∥W(B i , B i )∥ 2 F = k n C k n i W 2 ii + C k-2 n-2 i j,j̸ =i W 2 ij . Here, the set {B 1 , B 2 , • • • , B C k n } represents all possiblecombinations of the index vectors choosing k items from n without repetition.Lemma A.2. We have S + be the set of|S + | = b samples from [N ], drawn with replacement and uniformly at random. Then, ∀t, X t ∈ R n×n , we have:E ι t [∥ 1 b i∈S+ ∇f i (X t ) -∇f (X t )∥ 2 F ] = N -b b(N -1) E ι t [∥∇f i (X t ) -∇f (X t )∥ 2 F ].Proof. The proof is exactly the same as in Lemma 2.8 of [5] .Lemma A.3. The tangent space T X J of manifold constructed by X ⊤ JX = J, with X ∈ R n×n , is :EQUATIONwhere Y = t Ỹ with t is a positive scalar approaching 0.Proof. Assuming point X ∈ R n×n lies on manifold J , we have: h(X) = X ⊤ JX -J. Moving along Y ∈ R n×n in the tangent space of X, we obtain:h(X + Y) =(X + Y) ⊤ J(X + Y) -J =X ⊤ JX + X ⊤ JY + Y ⊤ JX + Y ⊤ JY -J ① =X ⊤ JY + Y ⊤ JX + Y ⊤ JY ② =tX ⊤ J Ỹ + t Ỹ⊤ JX + t 2 Ỹ⊤ J Ỹ where step ① uses X ⊤ JX = J; step ② uses Y = t Ỹ.Since t is a positive scalar approaching 0, we can ignore the higher-order term: t 2 Ỹ⊤ J Ỹ. According to the properties of the tangent space of any manifold, we have: h(X + Y) = 0, In other words, X ⊤ JY + Y ⊤ JX = 0, i.e. we obtain the defining equation for the tangent space:T X J ≜ {Y ∈ R n×n | X ⊤ JY + Y ⊤ JX = 0}.B.1 On the Global Optimal Solution for Problem (7) In Section 2.1, we have demonstrated how to use the breakpoint search method to obtain an optimal solution for the case of V = ( c s s c ) of Problem (7) . Since the structure of the other three casesV ∈ {( c -s -s c ), ( -c -s s c ), ( c -s s -c)} is exactly the same except for the coefficients of Problem (8), we will provide the corresponding coefficients in Problem (8): min c,s a c +b s +c c2 +d c s +e s2 , and omit the specific analysis process.Case (a). V = ( c -s -s c ):a = P 11 + P 22 , b = -P 12 -P 21 , c = 1 2 ( Q11 + Q41 + Q14 + Q44 ), d = -1 2 ( Q21 + Q31 + Q12 + Q42 + Q13 + Q43 + Q24 + Q34 ), and e = 1 2 ( Q22 + Q32 + Q23 + Q33 ). Case (b). V = ( -c -s s c ):a = -P 11 + P 22 , b = -P 12 + P 21 , c = 1 2 ( Q11 -Q41 -Q14 + Q44 ), d = 1 2 ( Q21 -Q31 + Q12 -Q42 -Q13 + Q43 -Q24 + Q34 ), and e = 1 2 ( Q22 -Q32 -Q23 + Q33 ). Case (c). V = ( c -s s -c ):a = P 11 -P 22 , b = -P 12 + P 21 , c = 1 2 ( Q11 -Q41 -Q14 + Q44 ), d = 1 2 (-Q21 + Q31 -Q12 + Q42 + Q13 -Q43 + Q24 -Q34 ), and e = 1 2 ( Q22 -Q32 -Q23 + Q33 ).C Proofs for Section 2C.1 Proof of Lemma 2.1Proof. Defining J BB = J(U B , U B ) , then we have:EQUATIONU ⊤ B JU B = J BB . Part (a). For any V ∈ R 2×2 and B ∈ {B i } C 2 n i=1, we have:[X + ] ⊤ JX + -X ⊤ JX ① = X ⊤ JU B (V -I 2 )U ⊤ B X + [U B (V -I 2 )U ⊤ B X] ⊤ JX +[U B (V -I 2 )U ⊤ B X] ⊤ J[U B (V -I 2 )U ⊤ B X] = X ⊤ [JU B (V -I 2 )U ⊤ B + U B (V -I 2 ) ⊤ U ⊤ B J + U B (V -I 2 ) ⊤ U ⊤ B JU B (V -I 2 )U ⊤ B ]X = X ⊤ [U B J BB (V -I 2 )U ⊤ B + U B (V -I 2 ) ⊤ J BB U ⊤ B + U B (V -I 2 ) ⊤ J BB (V -I 2 )U ⊤ B ]X = X ⊤ U B [J BB (V -I 2 ) + (V -I 2 ) ⊤ J BB + (V -I 2 ) ⊤ J BB (V -I 2 )]U ⊤ B X = X ⊤ U B [V ⊤ J BB V -J BB ]U ⊤ B X ② = 0.Part (b). Using the update rule forX + = X + U B (V -I 2 )U ⊤ B X ∈ R n×n , we derive: ∥X + -X∥ F = ∥U B (V -I 2 )U ⊤ B X∥ F ① ≤ ∥U B ∥ F • ∥(V -I 2 )U ⊤ B X∥ F , ② ≤ ∥U B ∥ F • ∥(V -I 2 )∥ F • ∥U ⊤ B ∥ F • ∥X∥ F , ③ = ∥V -I 2 ∥ F • ∥X∥ F ,where step ① and step ② use the norm inequality that ∥AX∥ F ≤ ∥A∥ F • ∥X∥ F for any A and X;step ③ uses ∥U B ∥ = ∥U ⊤ B ∥ = 1.Part (c). We define Z ≜ U ⊤ B X. We derive:∥X + -X∥ 2 H = ∥U B (V -I 2 )Z∥ 2 H ① = vec(U B (V -I 2 )Z) ⊤ Hvec(U B (V -I 2 )Z) ② = vec(V -I 2 ) ⊤ (Z ⊤ ⊗ U B ) ⊤ H(Z ⊤ ⊗ U B )vec(V -I 2 ) = ∥V -I 2 ∥ 2 (Z ⊤ ⊗UB) ⊤ H(Z ⊤ ⊗UB) ③ ≤ ∥V -I 2 ∥ 2 Q , where step ① uses ∥X∥ 2 H = vec(X) ⊤ Hvec(X); step ② uses (Z ⊤ ⊗ R)vec(U) = vec(RUZ) for all R, Z and U of suitable dimensions; step ③ uses the choice of Q ≽ Q ≜ (Z ⊤ ⊗ U B ) ⊤ H(Z ⊤ ⊗ U B ). C.2 Proof of Lemma 2.3Proof. We denote w = c + e. According to the properties of trigonometric functions, we have: (i)c2 = 1 1- t2 ; (ii) s2 = t2 1- t2 ; (iii) t = s c , leading to: c = ±1 √ 1- t2 , s = ± t √ 1- t2 with | t | < 1.We discuss two cases for Problem (8) .Case (a). c = 1 √ 1- t2 , s = t √ 1- t2 . Problem (8) is equivalent to the following problem: μ+ = arg min µ a+ t b √ 1- t2 + w+ t d 1-t2 -e. Therefore, the optimal solution μ+ can be computed as:cosh(μ + ) = 1 √ 1-( t+) 2 , and sinh(μ + ) = t+ √ 1-( t+) 2 (14) Case (b). c = -1 √ 1- t2 , s = -t √ 1- t2 . Problem (8) is equivalent to the following problem: μ-= arg min µ -a-t b √ 1- t2 + w+ t d 1-t2 -e. Therefore, the optimal solution μcan be computed as:EQUATION)We define the objective function as: F (c, s) ≜ ac + bs + cc 2 + dcs + es 2 . In view of ( 14) and ( 15), the optimal solution pair [cosh(μ, sinh(μ)] for problem (8) can be computed as:[cosh(μ), sinh(μ)] = arg min [c,s] F (c, s), s. t. [c, s] ∈ {[cosh(μ + ), sinh(μ + )], [cosh(μ -), sinh(μ -)]}Importantly, it is not necessary to compute the values μ+ for ( 14) and μfor (15) .C.3 Proof of Lemma 2.4Proof. The objective function for B t (i) as in Equation ( 3) is formulated as :f (X t ) + 1 2 ∥V i -I∥ 2 Q+θI + ⟨V i -I, [∇f (X t )(X t ) T ] B t (i) B t (i) ⟩ Part (1). For the part of 1 2 ∥V i -I∥ 2 Q+θI , it is obviously irrelevant. Part (2). For the part of ⟨V i -I, [∇f (X t )(X t ) T ] B t (i) B t (i) ⟩, we note that [∇f (X t )(X t ) ⊤ ] B t (i) B t (i) = [∇f (X t )](B t (i) , :)[(X t ) ⊤ ](:, B t (i) ) = [∇f (X t )](B t (i) , :)[(X t )(B t (i) , :)] ⊤ , which just use the informa- tion of block B t (i). The proof ends.C.4 Proof of Lemma 2.5Proof. Part (a). For the purpose of analysis, we define the following:∀i ∈ [ n 2 ], K i = U B (i) (V i - I 2 )U ⊤ B (i) X. ∥ n 2 i=1 [U B (i) (V i -I 2 )U ⊤ B (i) X]∥ 2 F ① =     K 1 K 2 . . . K n 2     2 F ② = ∥K 1 ∥ 2 F + ∥K 2 ∥ 2 F + • • • + ∥K n 2 ∥ 2 F ③ = n 2 i=1 [∥U B (i) (V i -I 2 )U ⊤ B (i) X∥ 2 F ]where step ① uses the definition of K i and the assumption that B ∈ Υ; step ② uses the definition of Squared Frobenius Norm; step ③ uses the definition of K i .Part (b). Using the update rule forX + = X + [ n/2 i=1 U B (i) (V i -I 2 )U ⊤ B (i) ]X ∈ R n×n ,we have the following inequalities:EQUATIONEQUATIONEQUATIONwhere step ① uses the conclusion of Part (a); step ② uses the same proof process of Part (b) of lemma 2.1.Part (c). We derive the following results:1 2 ∥X + -X∥ 2 H = 1 2 ∥[ n/2 i=1 U B (i) (V i -I 2 )U ⊤ B (i) ]X∥ 2 H ① = 1 2 n/2 i=1 ∥[U B (i) (V i -I 2 )U ⊤ B (i) ]X∥ 2 H ② ≤ 1 2 n/2 i=1 ∥V i -I 2 ∥ 2 Qwhere step ① uses the conclusion of Part (a); step ② uses the same proof process of Part (c) of lemma 2.1.Part (d). We derive the following results:EQUATION)where step ① uses ∀A, B, Proof. We consider the Lagrangian function of problem (1):1 2 ∥A -B∥ 2 F = 1 2 ∥A∥ 2 F + 1 2 ∥B∥ 2 F -⟨A, B⟩ ≥ 0, with A = ∥X + -X∥ 2 F and B = ∥[∇f (X) -G]∥ 2 F ; step ②EQUATION)Setting the gradient of L(X, Λ) w.r.t. X to zero yields:EQUATIONPart (a). Multiplying both sides by X T and using the fact that X T JX = J, we have JΛ = X T ∇f (X). Multiplying both sides by J T and using J ⊤ J = I, we have Λ = JX ⊤ ∇f (X). Since Λ is symmetric, we have Λ = ∇f (X) T XJ. Putting this equality into Equality (21) yields the following first-order optimality condition for Problem (1):EQUATIONPart (b). We let G = ∇f (X). We derive the following results:G = JXG T XJ ① ⇒ JX T • G = JX T • JXG T XJ ② ⇒ JX T G = G T XJ ③ ⇒ X(JX T G)X T = X(G T XJ)X T ④ ⇒ X JX T GX T J ≜G T J = J JXG T XJ ≜G X T (23) ⑤ ⇒ (XG T J) • JX = (JGX T ) • JX ⑥ ⇒ XG T X = JGJ ⑦ ⇒ JXG T XJ = G,where step ① uses the results of left-multiplying both sides by JX T ; step ② uses J•X T JX = JJ = I;step ③ uses the results of left-multiplying both sides by X and subsequently right-multiplying them by X T ; ④ uses G = JXG T XJ; step ⑤ uses the the results of right-multiplying both sides by JX;step ⑥ uses JJ = I and X T JX = J; step ⑦ uses the results of left-multiply both sides by J and right-multiplied by J.Given Equality (23) , we conclude that the critical point condition is equivalent to the requirement that the matrix X∇f ( X) T J is symmetric, which is expressed asXG T J = [XG T J] T .Proof. We use Ẍ and X to denote any BS-point and critical point, respectively.For all B ∈ Ω ≜ {B 1 , B 2 , . . . , B C 2 n }, we have:I 2 ∈ arg min V∈JB G(V; Ẍ, B).whereG(V; X, B) ≜ f (X) + 1 2 ∥V -I 2 ∥ 2 Q+θI + ⟨V -I, [∇f (X)(X) T ] BB ⟩.The Euclidean gradient of G(V; Ẍ, B) can be computed as:EQUATIONGiven Lemma 3.1, we set the Riemannian gradient of G(V; Ẍ, B) w.r.t. V to zero, leading to the following first-order optimality condition:EQUATIONLetting V = I 2 , and using the definition of G, we have:0 2,2 = [∇f (X)(X) T ] BB -J BB G⊤ J BB , ∀B ∈ {B i } C 2 n i=1 ⇒ 0 2,2 = U ⊤ B [∇f ( Ẍ) ẌT ]U B -J BB U ⊤ B [ Ẍ∇f ( Ẍ) T ]U B J BB , ∀B ∈ {B i } C 2 n i=1 ① ⇒ 0 2,2 = U ⊤ B [∇f ( Ẍ) ẌT ]U B -U ⊤ B J[ Ẍ∇f ( Ẍ) T ]JU B , ∀B ∈ {B i } C 2 n i=1 ② ⇒ 0 n,n = [∇f ( Ẍ) ẌT ] -J[ Ẍ∇f ( Ẍ) T ]J, ③ ⇒ [J∇f ( Ẍ) ẌT ] = [J∇f ( Ẍ) ẌT ] T ,where step ① uses U ⊤ B J = J BB U ⊤ B and JU B = U B J BB ; step ② uses the the following results for anyEQUATIONstep ③ uses the fact that both sides are left-multiplied by J. We conclude that the matrix J∇f ( Ẍ) ẌT is symmetric. Using Claim (b) of Lemma 3.1, we conclude that Ẍ is a also a critical point.Notably, the condition in Equation ( 25) is a necessary but not sufficient condition. This is because BS-point is the global minimum of Problem: arg min V∈JB G(V; Ẍ, B), according to Definition 3.2.E.1 Proof of Lemma 4.5Proof. By the definition of Gt , we haveE ι t [∥ Gt -∇f (X t )∥ 2 F ] ① = pE ι t [∥ 1 b b i=1 ∇f i (X t ) -∇f (X t )∥ 2 F ] + (1 -p)E ι t [∥ Gt-1 + 1 b ′ b ′ i=1 (∇f i (X t ) -∇f i (X t-1 )) -∇f (X t )∥ 2 F ] ② = pE ι t [∥ 1 b b i=1 ∇f i (X t ) -∇f (X t )∥ 2 F ] + (1 -p)E ι t-1 [∥ Gt-1 -∇f (X t-1 )∥ 2 F ] +(1 -p)E ι t [∥ 1 b ′ b ′ i=1 (∇f i (X t ) -∇f i (X t-1 )) -∇f (X t ) + ∇f (X t-1 )∥ 2 F ]where step ① uses formula (9); step ② uses that Gt-1 -∇f (X t-1 ) is measurable w.r.t. ι t-1 andEQUATIONwhere step ① uses that for any random variable X, Proof. For simplicity, we use B instead of B t . We will show that the following inequality holds :E[(X -E[X]) 2 ] ≤ E[X 2 ]; step ②EQUATION)Since Vt is the global optimal solution of Problem ( 5), we have:G( Vt ; X t , B) ≤ G(V; X t , B), V ∈ J B Letting V = I 2 , we have: G( Vt ; X t , B) ≤ G(I 2 ; X t , B).We further obtain:EQUATION)Using Inequality (2) with N = 1 and Part (c) of Lemma 2.1, we have:EQUATION)Adding Inequality ( 29) and ( 30) together, we obtain the inequality in (28) . Using the result of Part (b) in Lemma 2.1 that∥X + -X∥ 2 F ∥X∥ 2 F ≤ ∥V -I 2 ∥ 2F , we have the following sufficient decrease condition:EQUATION)We now prove the global convergence. Taking the expectation for Inequality (31) , we obtain a lower bound on the expected progress made by each iteration for Algorithm 1:E ξ t+1 [f (X t+1 )] -E ξ t [f (X t )] ≤ -E ξ t [ θ 2 ∥ Vt -I 2 ∥ 2 F ].Summing up the inequality above over t = 0, 1, . . . , T , we have:E ξ t [ θ 2 T t=0 ∥ Vt -I 2 ∥ 2 F ] ≤ f (X 0 ) -E ξ T +1 [f (X T +1 )] ≤ f (X 0 ) -f ( X).As a result, there exists an index t with 0 ≤ t ≤ T such thatEQUATIONFurthermore, for any t, we have:EQUATIONCombining Inequality (32) and equality (33) , we have the following result:EQUATION)We will give the arithmetic operations of GS-JOBCD. By the chosen parameters and Inequality (34), we haveE(X t) ≤ 2(f (X 0 )-f ( X)) θ(T +1) ≤ ϵ.We define ∆ 0 = f (X 0 ) -f ( X) and set T + 1 = 2∆0 ϵθ . Denoting m t to be the number of arithmetic operations at t-th iteration, we have for t ≥ 1:E ξ t [m t ] = O(2N ).Then we have for t ≥ 1, the total number of arithmetic operations M T in T iterations to obtain ϵ-BS-point isE ξ T [M T ] = E ξ t [ T t=0 m t ] = 2(T + 1)N = O((T + 1)N ). We have (T + 1)N = N 2∆0 ϵθ = O( ∆0N ϵ ).E.3 Proof of Theorem 4.7Proof. For simplicity, we use B instead of B t . Defining V: t as the global optimal solution of arg min V: T (V : ; X t , B), we have:T ( V: t ; X t , B) ≤ T (V : ; X t , B), ∀i, V i ∈ J B (i)Letting V i = I 2 , ∀i, we have: T ( V: t ; X t , B) ≤ T (I 2 ; X t , B). We further obtain:EQUATION)Using the results of telescoping Inequality (2) over i from 1 to N with Part (c) of Lemma 2.5, we have:EQUATIONAdding inequality (35) , and (36) together, we obtain the inequality in (37) .EQUATION)where step ① uses Part (d) of Lemma 2.5.Taking expectation on both sides of inequality (37) with respect to all randomness of the algorithm, and adding the inequality in Lemma 4.5 × 1 2p to (37), we have:EQUATIONSumming up the inequality above over t = 0, 1, . . . , T , we have:EQUATIONAs a result, there exists an index t with 0 ≤ t ≤ T such thatEQUATIONDefiningϖ = θ-X 2 2 - L 2 f X 2 (1-p) 2pb ′, furthermore, for any t and ∀i, we have:EQUATIONCombining inequality ( 40) and ( 41) , we have the following result:EQUATIONBy the chosen parameters and Inequality (42), we haveE(X t) ≤ 1 (T +1)ϖ (f (X 0 ) -f ( X) + (T +1)(N -b) 2b(N -1) σ 2 + 1-p 2p (u 0 -E ι T +1 [u T +1 ])) ≤ ϵ.We define ∆ 0 = f (X 0 ) -f ( X) and set T + 1 = ∆0 ϵϖ . Denoting m i t to be the number of arithmetic operations to update the i-th block at t-th iteration, we have for t ≥ 1E ι t [m i t ] = O(2(pb + (1 -p)b ′ )). Letting m t be the number of arithmetic operations in the t-the iteration, we have for t ≥ 1E ι t [m t ] = E ι t [ n/2 i=1 m i t ] = O((pb + (1 -p)b ′ )n/2 × 2) = O(n(pb + (1 -p)b ′ )).Hence, the total number of arithmetic operations M T in T iterations to obtain ϵ-BS-point isE ι T [M ] = E ι t [ T t=0 m t ] = O(bn) + E ι t [ T t=1 m t ] = O(bn + T n(pb + (1 -p)b ′ )). Since b = N, b ′ = √ b and p = b ′ b+b ′ , ϖ = θ-X 2 2 - L 2 f X 2 (1-p) 2pb ′ = 1 2 (θ -X 2 -L 2 f X 2 ), we have nT (pb + (1 -p)b ′ ) = n ∆0 ϵ(θ-X 2 -L 2 f X 2 ) 2bb ′ b+b ′ ≤ n∆0 ϵ(θ-X 2 -L 2 f X 2 ) 2b ′ = O( ∆0 √ N ϵ).Proof. For simplicity, we use B instead of B t . We notice that the Riemannian gradient of T (V : ; X t , B)at the point V i = I 2 , ∀i . Defining G = Gt [X t ] ⊤ and using JU B = U B J BB , U ⊤ B J = J BB U ⊤ B ,wehave:EQUATIONThen, we prove the following important lemmas.Lemma E.1. We have the following result for VR-J-JOBCD:E ι t+1 [∥ Gt -Gt+1 ∥ F ] ≤ pE ι t [ √ u t ] + L f E ι t+1 [∥X t -X t+1 ∥ F ]Proof. By the definition of Gt , with the choice of b = N , b ′ = √ b and p = b ′ b+b ′ , we haveE ι t+1 [∥ Gt -Gt+1 ∥ F ] ① = E ι t+1 [∥ Gt -p b b i=1 ∇f i (X t+1 ) -1-p b ′ b ′ i=1 (∇f i (X t+1 ) -∇f i (X t )) -(1 -p) Gt ∥ F ] = E ι t+1 [∥p Gt -p b b i=1 ∇f i (X t+1 ) -1-p b ′ b ′ i=1 (∇f i (X t+1 ) -∇f i (X t ))∥ F ] ② ≤ pE ι t+1 [∥ Gt -∇f (X t+1 )∥ F ] + 1-p b ′ E ι t [∥ b ′ i=1 ∇f i (X t+1 ) -∇f i (X t )∥ F ] ③ ≤ pE ι t [∥ Gt -∇f (X t )∥ F ] + pE ι t+1 [∥∇f (X t ) -∇f (X t+1 )∥ F ] + 1-p b ′ E ι t+1 [∥ b ′ i=1 ∇f i (X t+1 ) -∇f i (X t )∥ F ] ④ ≤ pE ι t [ √ u t ] + pE ι t+1 [∥∇f (X t ) -∇f (X t+1 )∥ F ] + (1 -p)E ι t+1 [∥∇f i (X t+1 ) -∇f i (X t )∥ F ] ⑤ ≤ pE ι t [ √ u t ] + L f E ι t+1 [∥X t -X t+1 ∥ F ]where step ① uses formula (9); step ② uses norm inequality and1 b b i=1 ∇f i (X t+1 ) = ∇f (X t+1 )with b = N and norm inequality; step ③ uses triangle inequality that ∥A -B∥ F ≤ ∥A -C∥ F + ∥C -B∥ F , for any A, B and C; step ④ the definition of u t ; step ⑤ uses Inequality (2) and the results of telescoping it over i from 1 to N .Lemma E.2. (Riemannian gradient Lower Bound for the Iterates Gap) We de-fine ϕ ≜ (3X + VX)G + (1 + V 2 + n 2 (X 2 + V 2 X 2 ))L f + (1 + V 2 )θ.It holds that:E ι t+1 [dist(0, ∇ J T (I 2 ; X t+1 , B t+1 ))] ≤ ϕ • E ι t [ n/2 i=1 ∥ Vt i -I 2 ∥ F ] + np √ u t 2 (X + V 2 X).Proof. For notation simplicity, we define:EQUATIONEQUATIONEQUATIONFirst, using the optimality of Vt i , i ∈ {1, • • • , n 2 } for the subproblem, we have:EQUATIONwhereEQUATION)Using the relation that Gi = Υ i1 + Υ i2 , we obtain the following results from the above equality:EQUATIONwhere step ① uses Υ i2 = Ω i1 + Ω i2 . Then we derive the following results:EQUATIONwhere step ① uses Equality (43) ; step ② uses the fact that both the working set B t and B t+1 are selected randomly and uniformly; step ③ uses the definition of Ω i0 in (44); step ④ uses -Ω i1 + Ω i1 = 0 and -Ω ⊤ i1 + Ω ⊤ i1 = 0; step ⑤ uses the norm inequality; step ⑥ uses the norm inequality; step ⑦ uses the norm inequality; step ⑧ uses Equality (49); step ⑨ uses the norm inequality. We now establish individual bounds for each term for Inequality (50) .For the first termEQUATIONwhere step ① uses[X t -X t ] BiBi = U B (i) ( Vt i -I 2 )U ⊤ B (i) X t ; step ② uses the inequality ∥XY∥ F ≤∥X∥ F ∥Y∥ F for all X and Y repeatedly and the fact that ∀t, ∥ Gt ∥ F ≤ G and ∀t, ∥X t ∥ F ≤ X; step ③ uses the norm inequality.For the second term 50):E ι t [∥ n/2 i=1 J B (i) Vt i Υ ⊤ i1 Vt i J B (i) -Υ i1 ∥ F ] in (EQUATIONwhere step ① uses the triangle inequality; step ② uses the inequality ∥XY∥ F ≤ ∥X∥ F ∥Y∥ F for all X and Y and ∀t, ∥V t ∥ F ≤ V; step ③ uses the definition of Υ i1 ; step ④ uses the choice of Q ⪯ L f I and the norm inequality.For the third term (50) , we have:E ι t [∥ n/2 i=1 Vt i Ω ⊤ i1 Vt i -Ω ⊤ i1 ∥ F ] inEQUATIONwhere step ① uses the fact that-Vt i Ω ⊤ i1 I 2 + Vt i Ω ⊤ i1 = 0; step ② uses the norm inequality and ∀t, ∥V t ∥ F ≤ V; step ③ uses the fact that ∥Ω i1 ∥ F = ∥U ⊤ B (i) Gt [X t ] ⊤ U B (i) ∥ F ≤ X∥ Gt ∥ F , ∀i whichcan be derived using the norm inequality ; step ④ uses the fact that ∀X, ∥ Gt ∥ F ≤ G.For the fourth term 50), we have:E ι t [∥ n/2 i=1 J B (i) Vt i Ω ⊤ i2 Vt i J B (i) -Ω i2 ∥ F ] in (EQUATIONwhere step ① uses the triangle inequality; step ② uses the norm inequality and ∀t, ∥V t ∥ F ≤ V; step ③ uses the definition of ∀i, (46) ; step ④ uses the norm inequality and ∀t, ∥X t ∥ F ≤ X; step ⑤ uses Lemma E.1; step ⑥ uses Part (b) in Lemma 2.5 and ∀t, ∥X t ∥ F ≤ X.Ω i2 = U ⊤ B (i) [ Gt -Gt ][X t ] ⊤ U B (i) inIn view of( 51), ( 52), ( 53), (54) , and (50), we have:E ι t+1 [∥∇ J T (I 2 ; X t+1 , B t+1 )∥ F ] ≤ np 2 (X + V 2 X)E ι t [ √ u t ] + (c 1 + c 2 + c3 + c4) • E ι t [ n/2 i=1 ∥ Vt i -I 2 ∥ F ] = np 2 (X + V 2 X)E ι t [ √ u t ] + ϕE ι t [ n/2 i=1 ∥ Vt i -I 2 ∥ F ]wherec 1 = 2XG, c 2 = (1 + V 2 )(L f + θ), c 3 = (X + VX)G, and c 4 = n 2 (X 2 + V 2 X 2 )L f .Lemma E.3. We have the following results: dist(0,∇ J f (X t )) ≤ γ • ∥∇ J T (I 2 ; X t , B)∥ F + 2X 2 E ι t [u t ] with γ ≜ X C 2 n .Proof. We have the following inequalities:∥∇ J f (X t )∥ F ① = ∥∇f (X t ) -JX t (∇f (X t )) ⊤ X t J∥ F ② = ∥∇f (X t )(X t ) ⊤ JX t J -JX t (∇f (X t )) ⊤ JJX t J∥ F ③ ≤ ∥∇f (X t )(X t ) ⊤ -JX t (∇f (X t )) ⊤ J∥ F ∥JX t J∥ F ④ ≤ X∥∇f (X t )(X t ) ⊤ -JX t (∇f (X t )) ⊤ J∥ Fwhere step ① uses the definition of ∇ J f (X t ); step ② uses JJ = I and X ⊤ JX = J ⇒ X ⊤ JXJ = JJ = I; step ③ uses the norm inequality and ; step ④ uses ∀t, ∥X t ∥ F ≤ X.We Consider ∥∇f (X t )(X t ) ⊤ -JX t (∇f (X t )) ⊤ J∥ F : ∥∇f (X t )(X t ) ⊤ -JX t (∇f (X t )) ⊤ J∥ F ① ≤ ∥ Gt (X t ) ⊤ -JX t ( Gt ) ⊤ J∥ F + ∥(∇f (X t ) -Gt )(X t ) ⊤ -JX t (∇f (X t ) -Gt ) ⊤ J∥ F ② ≤ ∥ Gt (X t ) ⊤ -JX t ( Gt ) ⊤ J∥ F + ∥∇f (X t ) -Gt ∥ F • ∥X t ∥ F + ∥X t ∥ F • ∥∇f (X t ) -Gt ∥ F ③ ≤ ∥ Gt (X t ) ⊤ -JX t ( Gt ) ⊤ J∥ F + 2X E ι t [u t ]where step ① uses ∀A, B, ∥A∥ F -∥B∥ F ≤ ∥A -B∥ F ; step ② uses the norm inequality; step ③ uses ∀t, ∥X t ∥ F ≤ X. Thus,∥∇ J F (X t )∥ F ≤ X∥ Gt (X t ) ⊤ -JX t ( Gt ) ⊤ J∥ F + 2X 2 E ι t [u t ] ① ≤ X C 2 n • ∥ n/2 i=1 U ⊤ B (i) [ Gt (X t ) ⊤ -JX t ( Gt ) ⊤ JU B (i) ∥ F ] + 2X 2 E ι t [u t ] ② = X C 2 n • ∥∇ J T (I 2 ; X t , B)∥ F + 2X 2 E ι t [u t ]where step ① uses Lemma A.1 with W = Gt (X t ) ⊤ -JX t ( Gt ) ⊤ J and k = 2; step ② uses the definition of ∇ J T (I 2 ; X t , B).We now present the following useful lemma.Lemma E.4. We defineT X J ≜ {Y ∈ R n×n | A X (Y) = 0} and A X (Y) ≜ X ⊤ JY + Y ⊤ JX.For any G ∈ R n×n and X ⊤ JX = J, the unique minimizer of the following optimization problem:Ȳ = arg min Y∈T XM J h(Y) = 1 2 ∥Y -G∥ 2 F , satisify h( Ȳ) ≤ h(G -JXG ⊤ XJ).Proof. We note that Ȳ = arg min Y∈T X J1 2 ∥Y -G∥ 2 F = arg min Y 1 2 ∥Y -G∥ 2 F , s.t. X ⊤ JY + Y ⊤ JX = 0.Introducing a multiplier Λ ∈ R n×n for the linear constraints X ⊤ JY + Y ⊤ JX = 0, we have following Lagrangian function:L(Y; Λ) = 1 2 ∥Y -G∥ 2 F + ⟨X ⊤ JY + Y ⊤ JX, Λ⟩.We naturally derive the following first-order optimality condition: Y -G + JXΛ = 0, X ⊤ JY + Y ⊤ JX = 0. Incorporating the term Y = G -JXΛ into X ⊤ JY + Y ⊤ JX = 0, we obtain:EQUATIONAny Λ satisfying formula ( 55) is a feasible point, so we can easily find :EQUATION)where step ① uses the fact that any matrix X satisfying the J-orthogonality constraint has a determinant of 1 or -1, thus inv(X) exists; step ② multiply both sides of the equation by XJ;step ③ uses X T JX = J and JJ = I; step ④ multiply both sides of the equation by J and uses JJ = I; step ⑤ uses the fact that Λ is a symmetric matrix.Therefore, a feasible solution Y can be computed asY = G -JXΛ = G -JXG ⊤ XJ. Since Ȳis the optimal solution, there must be h( Ȳ) ≤ h(G -JXG ⊤ XJ).We now present the proof of this lemma.Lemma E.5. For any X ∈ R n×n , it holds that dist(0, ∇f • (X)) ≤ dist(0, ∇ J f (X)).Proof. For the purpose of analysis, we define the nearest J orthogonal matrix to an arbitrary matrix Y ∈ R n×n is given by P J (X). Similarly, we have P T X J (∇f (X)) for projecting gradient ∇f (X) into space T X J .We recall that the following first-order optimality conditions are equivalent for all X ∈ R n×n :EQUATIONTherefore, we derive the following results:EQUATIONWe let G ∈ ∇f (X) and obtain the following results from the above equality:EQUATIONEQUATION)where step ① uses Lemma E.4; step ② uses∇ J f (X) = G -JXG ⊤ XJ with G ∈ ∇f (X).First of all, since f • (X) ≜ f (X) + I J (X) is a KL function, we have from Proposition 4.8 that:EQUATIONwhere step ① uses Lemma E.5. Here, φ(•) is some certain concave desingularization function. Since φ(•) is concave, we have:EQUATIONApplying the inequality above with ∆ = f (X t ) -f ( X) and ∆ + = f (X t+1 ) -f ( X), we have:EQUATION)With the sufficient descent condition as shown in Theorem 4.7, we derive the following inequalities:EQUATIONwhere step ① uses ∀t, ∥X t ∥ F ≤ X.EQUATIONwhere step ① uses the sufficient descent condition as shown in Theorem 4.7; step ② uses Inequality (64) and (62) with X ′ = X t and X = X; step ③ uses lemma E.3 ; step ④ uses Lemma E.2 ; step ⑤ uses∀x i ∈ R, x1+•••+xn n ⩽ x 2 1 +•••+x 2 n n ; step ⑥ applies the inequality that ∀θ ′ > 0, a, b, ab ≤ θ ′ a 2 2 + b 2 2θ ′ with a = n/2 i=1 ∥ Vt-1 i -I 2 ∥ 2 F , b = E t γϕ n 2 ; a = E ι t [u t ], b = E t (2X 2 + γ np 2 X + γ np 2 V 2 X); step ⑦ denote A 2 ≜ (2X 2 +γ np 2 X+γ np 2 V 2 X) 2 2 θ+ nγ 2 ϕ 2 4θ ′ . To simplify the formula, we defineℵ t = n/2 i=1 ∥ Vt i -I 2 ∥ 2 F .Multiplying both sides by 2 and taking the square root of both sides, we have:EQUATIONTo recursively eliminate term ( θ + 1)E ι t [u t ], we take the root of both sides of the Inequality in Lemma 4.5:E ι t [u t ] ≤ p(N -b) b(N -1) σ 2 + (1 -p)E ι t-1 [u t-1 ] + L 2 f X 2 (1-p) b ′ E ι t-1 [ℵ t-1 ] ≤ p(N -b) b(N -1) σ 2 + (1 -p) E ι t-1 [u t-1 ] + L 2 f X 2 (1-p) b ′ E ι t-1 [ℵ t-1 ] (70) Adding Inequality √ θ+1 1- √ 1-p × (70) to (69) E ι t [ θ -X 2 √ ℵ t ] ≤ E t A + ( √ θ ′ + L 2 f X 2 (1-p) b ′ √ θ+1 1- √ 1-p )E ι t-1 [ √ ℵ t-1 ] + √ 1-p √ ( θ+1) 1- √ 1-p ( E ι t-1 [u t-1 ] -E ι t [u t ]) + √ θ+1 1- √ 1-p p(N -b) b(N -1) σ 2 (71)With the choice√ θ ′ = √ θ-X 2 2 - L 2 f X 2 (1-p) b ′ √ θ+1 1- √1-p , we have:EQUATION2)E ι t-1 [ √ ℵ t-1 ] + √ 1-p √ ( θ+1) 1- √ 1-p ( E ι t-1 [u t-1 ] -E ι t [u t ]) + √ θ+1 1- √ 1-p p(N -b) b(N -1) σ 2 (72)Rearranging terms, we have:EQUATIONSumming the inequality above over t = 1, 2 . . . , T , we have:E ι T [ θ -X 2 √ ℵ T ] + E ι T -1 [ √ θ-X 2 2 T -1 t=1 √ ℵ t ] ≤ A T t=1 E t + √ 1-p √ ( θ+1) 1- √ 1-p ( E ι 0 [u 0 ] -E ι T [u T ]) + T √ θ+1 1- √ 1-p p(N -b) b(N -1) σ 2 + √ θ-X 2 2 √ ℵ 0 ≤ A T t=1 E t + √ 1-p √ ( θ+1) 1- √ 1-p N -b b(N -1) σ 2 + T √ θ+1 1- √ 1-p p(N -b) b(N -1) σ 2 + E ι t [ √ θ-X 2 2 √ ℵ 0 ] ≤ A T t=1 E t + √ 1-p √ ( θ+1) 1- √ 1-p N -b b(N -1) σ 2 + T √ θ+1 1- √ 1-p p(N -b) b(N -1) σ 2 + √ θ-X 2 2 n 2 (V + √ 2) 2where step ① uses the fact thatE ι T [u T ] ≥ 0 and E ι 0 [u 0 ] ≤ N -b b(N -1) σ 2 ; step ② uses ∀t, ∥V∥ F ≤ V, then, ∥V i -I 2 ∥ 2 F ≤ (∥V i ∥ F + ∥I 2 ∥ F ) 2 ≤ (X + √ 2) 2 and n/2 i=1 ∥ V0 i -I 2 ∥ 2 F ≤ n 2 (V + √ 2) 2 . Define C = √ 1-p √ ( θ+1) 1- √ 1-p N -b b(N -1) σ 2 + T √ θ+1 1- √ 1-p p(N -b)b(N -1) σ 2 and rearrange terms, we have:EQUATIONT t=1 E t , we have:EQUATIONwhere step ① uses the definition of E i in (64); step ② uses a basic recursive reduction; step ③ uses the fact the desingularization function φ(•) is positive. Combining Inequality (74) and ( 75), we obtain :E ι t [ θ-X 2 2 T t=1 √ ℵ t ] ≤ Aφ(f (X 1 ) -f ( X)) + C + √ θ-X 2 2 n 2 (V + √ 2) 2Using the inequality that∥X + -X∥ 2 F X 2 ≤ ∥X + -X∥ 2 F ∥X∥ 2 F ≤ n/2 i=1 ∥ Vi -I 2 ∥ 2 F as shown in Part (b) in Lemma 2.5, we have:E ι t [ θ-X 2 2X T t=1 ∥X t+1 -X t ∥ F ] ≤ Aφ(f (X 1 ) -f ( X)) + C + √ θ-X 2 2 n 2 (V + √ 2) 2 Since b = N, b ′ = √ b and p = b ′ b+b ′ , C = √ 1-p √ ( θ+1) 1- √ 1-p N -b b(N -1) σ 2 + T √ θ+1 1- √ 1-p p(N -b)b(N -1) σ 2 = 0, we can get the expression for C:E ι t [ t j=1 ∥X j+1 -X j ∥ F ] ≤ C C ≜ 2X θ-X 2 (Aφ(f (X 1 ) -f ( X)) + √ θ-X 2 2 n 2 (V + √ 2) 2 )Considering that: 4 ), we have:EQUATIONA = (2X 2 +γ np 2 X+γ np 2 V 2 X) 2 2 θ + nγ 2 ϕ 2 4θ ′ = O( 1 N 1/4). Finally, we haveC = O( φ(f (X 1 )-f ( X)) N 1/4) E.5 Proof of Theorem 4.9Proof. For simplicity, we use B instead of B t . Initially, we prove the following important lemmas.Lemma E.6. (Riemannian gradient Lower Bound for the Iterates Gap) We define ϕ≜ (3X + VX)G + (1 + X 2 + V 2 + V 2 X 2 )L f + (1 + V 2 )θ. It holds that: E ξ t+1 [dist(0, ∇ J G(I 2 ; X t+1 , B t+1 ))] ≤ ϕ • E ξ t [∥ Vt -I 2 ∥ F ].Proof. The proof process is exactly the same as in lemma E.2 and will not be repeated here.The following lemma is useful to outline the relation of ∥∇ J f (X t )∥ F and ∥∇ J G(I 2 ; X t , B)∥ F .Lemma E.7. We have the following results:dist(0, ∇ J f (X t )) ≤ γ • E ξ t-1 [dist(0, ∇ J G(I 2 ; X t , B))] with γ ≜ X C 2 n .Proof. We have the following inequalities:∥∇ J f (X t )∥ 2 F ① = ∥G t -JX t (G t ) ⊤ X t J∥ 2 F ② = ∥G t (X t ) ⊤ JX t J -JX t (G t ) ⊤ JJX t J∥ 2 F ③ ≤ ∥G t (X t ) ⊤ -JX t (G t ) ⊤ J∥ 2 F ∥JX t J∥ 2 F ④ ≤ X 2 ∥W∥ 2 F , with W ≜ G t (X t ) ⊤ -JX t (G t ) ⊤ J ⑤ ≤ X 2 C 2 n • E ξ t-1 [∥U ⊤ B [G t (X t ) ⊤ -JX t (G t ) ⊤ J]U B ∥ 2 F ] ⑥ = X 2 C 2 n • E ξ t-1 [∥∇ J G(I 2 ; X t , B)∥ 2 F ]where step ① uses the definition of ∇ J f (X t ); step ② uses JJ = I and X ⊤ JX = J ⇒ X ⊤ JXJ = JJ = I; step ③ uses the norm inequality and ; step ④ uses the definition of W ≜ G t (X t ) ⊤ -JX t (G t ) ⊤ J and ∀t, ∥X t ∥ F ≤ X ; step ⑤ uses Lemma (A.1) with k = 2; step ⑥ uses the definition of ∇ J G(I 2 ; X t , B). Taking the square root of both sides, we finish the proof of this lemma.Finally, we obtain our main convergence results. First of all, since f • (X) ≜ f (X) + I J (X) is a KL function, we have from Proposition 4.8 that:EQUATIONwhere step ① uses Lemma E.5. Here, φ(•) is some certain concave desingularization function. Since φ(•) is concave, we have:∀∆ ∈ R, ∆ + ∈ R, φ(∆ + ) + (∆ -∆ + )φ ′ (∆) ≤ φ(∆).Applying the inequality above with ∆ = f (X t ) -f ( X) and ∆ + = f (X t+1 ) -f ( X), we have:EQUATION)We derive the following inequalities:E ξ t [ θ 2 ∥ Vt -I 2 ∥ 2 F ] ① ≤ E ξ t [f (X t ) -f (X t+1 )] ② ≤ E ξ t [ E t φ ′ (f (X t )-f ( X)) ] ③ ≤ E ξ t [E t ∥∇ J f (X t )∥ F ] ④ ≤ E ξ t [E t γ∥∇ J G(I 2 ; X t , B)∥ F ] ⑤ ≤ E ξ t-1 [E t γϕ| Vt-1 -I 2 ∥ F ] ⑥ ≤ E ξ t-1 [ θ ′ 2 ∥ Vt-1 -I 2 ∥ 2 F + (E t γϕ) 2 2θ ′ ],∀θ ′ > 0, where step ① uses the sufficient descent condition as shown in Theorem 4.6; step ② uses Inequality (77); step ③ uses Inequality (76) with X ′ = X t and X = X; step ④ uses Lemma E.7; step ⑤ uses Lemma E.6; step ⑥ applies the inequality that ∀θ ′ > 0, a, b, ab≤ θ ′ a 2 2 + b 2 2θ ′ with a = ∥ Vt-1 -I 2 ∥ F and b = E t γϕ.Multiplying both sides by 2 and taking the square root of both sides, we have:√ θE ξ t [∥ Vt -I 2 ∥ F ] ≤ (E t γϕ) 2 θ ′ + θ ′ E ξ t-1 [∥ Vt-1 -I 2 ∥ 2 F ], ∀θ ′ > 0 ① ≤ √ θ ′ E ξ t-1 [∥ Vt-1 -I 2 ∥ F ] + E t γϕ √ θ ′ , ∀θ ′ > 0,where step ① uses the inequality that √ a + b ≤ √ a + √ b for all a ≥ 0 and b ≥ 0. Summing the inequality above over i = 1, 2 . . . , t, we have:√ θE ξ t [∥ Vt -I 2 ∥ F ] - √ θ ′ E ξ 0 [∥ V0 -I 2 ∥ F ] + t-1 i=1 ( √ θ - √ θ ′ )E ξ i [∥ Vi -I 2 ∥ F ] ≤ γϕ √ θ ′ t i=1 E i ① = γϕ √ θ ′ t i=1 φ(f (X i ) -f ( X)) -φ(f (X i+1 ) -f ( X)) ② = γϕ √ θ ′ [φ(f (X 1 ) -f ( X)) -φ(f (X t+1 ) -f ( X))] ③ ≤ γϕ √ θ ′ φ(f (X 1 ) -f ( X)), where step ① uses the definition of E i in (77); step ② uses a basic recursive reduction; step ③ uses the fact the desingularization function φ(•) is positive. With the choice θ ′ = θ 4 , we have:EQUATIONEQUATIONwhere step ① uses ∀t, ∥V∥ F ≤ V,then, ∥V -I 2 ∥ F ≤ ∥V∥ F + ∥I∥ F ≤ V + √ 2.Finally, we obtain from Inequality (79):1 2 t i=1 E ξ i [∥ Vi -I 2 ∥ F ] ≤ 2γϕ θ φ(f (X 1 ) -f ( X)) + 1 2 (V + √ 2) ① ⇒ 1 2 t i=1 E ξ i [∥X i+1 -X i ∥ F ] ≤ ( 2Xγϕ θ φ(f (X 1 ) -f ( X)) + X 2 (V + √ 2))where step ① uses the inequality that∥X i+1 -X i ∥ F X ≤ ∥ Vi -I 2 ∥ F as shown in Part (b) in Lemma 2.1. Finally, we can get the expression for C:C ≜ 4Xγϕ θ φ(F (X 1 ) -F ( X)) + X(V + √ 2) = nO(φ(f (X 1 ) -f ( X)))To begin with, we give the definition of the Ultrahyperbolic manifold U p,q α , which will be used in Ultra-hyperbolic geodesic distance d α (x, y) and Diffeomorphism φ(•).▶ Ultrahyperbolic manifold. Vectors in an ultrahyperbolic manifold is defined asU p,q α = {x = (x 1 , x 2 , • • • , x p+q ) ⊤ ∈ R p,q : ∥x∥ 2 q = -α 2 }[48], where α is a non-negative real number denoting the radius of curvature. ∥x∥ 2 q = ⟨x, x⟩ q , ∀x, y ∈ R p,q , ⟨x, y⟩ q = p i=1 x i y i -p+q j=p+1 x j y j is a norm of the induced scalar product. The hyperbolic and spherical manifolds can be defined as:H α = U p,1 α , S α = U 0,q α .▶ Ultra-hyperbolic geodesic distance. The ultra-hyperbolic geodesic distance [27][28] d γ (•, •) is formulated: ∀x ∈ U p,q α , y ∈ U p,q α and α > 0, d α (x, y) = { α cosh -1 (| ⟨x,y⟩q α 2 |) if | ⟨x,y⟩q α 2 | ≥ 1 α cos -1 (| ⟨x,y⟩q α 2 |) otherwise. . ▶ Diffeomorphism.[Theorem 1 Diffeomorphism of [49] ]: Any vector x ∈ R p × R q * can be mapped into U p,q α by a double projection φ = ϕ -1 •ϕ, with ψ(x) = (s α t ∥t∥ ), ψ -1 (z) = ( v √ α 2 +∥v∥ 2 α u ), where x = ( s t ) ∈ U p,q α with s ∈ R p and t ∈ R q * • z = ( v u ) ∈ R p × S q α with v ∈ R p and u ∈ S q α .The J orthogonal matrix can be used as an isometric linear operator in the Ultrahyperbolic manifold, [48] et al. extended the knowledge graph model from hyperbolic space to Ultra-hyperbolic space (named as UltraE) by this property. The UltraE model is formulated as follows: The proposed methods generally give the best performance.min R,E,b L(R, E, b) ≜ -1 N (h,r,t)∈∆ (log s(h, r, t) + (h ′ ,r ′ ,r ′ )∈∆ ′ (h,r,t) log(1 -s(h ′ , r ′ , t ′ ))) s.t. s(h, r, t) = σ(-d 2 α (R r E h , E t ) + b h + b t + δ) R ⊤ r JR r = J where E ∈ R ne×n with E h = E(h, :) ∈ U p,q α , b ∈ R nr with b h = b(r) ∈ R, R ∈ R nr×n×n with R r = R(▶ Hyperbolic Structural Probe Problem. Table 3 and Figure J-JOBCD and VR-J-JOBCD often greatly improves upon CSDM, thus they show better MRR and hits results. (ii) In models with limited sample sizes, the computational efficiency of VR-J-JOBCD is inferior to that of J-JOBCD. This discrepancy arises because each iteration in VR-J-JOBCD necessitates two instances of backpropagation, thus consuming substantial computational resources.(iii) The proposed methods generally give the best performance.Table 2 : The convergence curve of the compared methods for solving HEVP. (+) indicates that after the convergence of the CSDM, UMCM and ADMM, utilizing the GS-JOBCD for optimization markedly enhances the objective value. The 1 st , 2 nd , and 3 rd best results are colored with red, green and blue, respectively. (n, p) represents the dimension and p-value of the J orthogonal matrix (square matrix). The value in () stands for (1000-100-90) -6.42e+03(1.2e-09) -6.42e+03(1.2e-09) -1.59e+04(1.2e-09) -4.51e+04(1.4e-08) -5.05e+04(5.6e-08) -2.64e+04(6.3e-09)(+) -3.35e+04(7.0e-09)(+) -3.89e+04(6.7e-09)(+) CnnCaltech (2000-1000-900) -3.10e+02(1.1e-08) -3.10e+02(1.2e-10) -3.41e+02(1.3e-10) -3.64e+02(1.7e-10) -1.65e+03(8.0e-09) -3.29e+02(1.1e-08)(+) -3.32e+02(1.3e-10)(+) -3.61e+02(1.4e-10)(+) gisette (3000-1000-900) -1.74e+06(1.2e-10) -1.74e+06(1.2e-10) -2.05e+06(1.2e-10) -2.57e+06(1.8e-10) -6.46e+06(8.0e-09) -2.00e+06(1.3e-10)(+) -1.99e+06(1.2e-10)(+) -2.33e+06(1.4e-10)(+) mnist (1000-780-650) -5.19e+04(1.9e-10) -5.19e+04(1.9e-10) -6.12e+04(2.2e-10) -1.02e+05(2.9e-10) -4.03e+05(1.2e-08) -6.75e+04(2.1e-10)(+) -6.58e+04(2.1e-10)(+) -n ij |X ⊤ JX -J| ij . dataname (m-n-p) UMCM ADMM CSDM GS-JOBCD J-JOBCD UMCM+GS-JOBCD ADMM+GS-JOBCD CSDM+GS-JOBCDEQUATION(100-100-90) -6.14e+03(1.1e-07) -6.14e+03(1.2e-09) -8.14e+03(1.4e-09) -8.31e+03(1.5e-08) -1.31e+04(5.8e-08) -7.74e+03(1.2e-07)(+) -7.77e+03(1.3e-08)(+) -9.69e+03(1.1e-08)(+) randn1000 (1000-1000-900) -6.33e+05(1.2e-10) -6.33e+05(1.2e-10) -6.84e+05(1.3e-10) -6.46e+05(1.9e-10) -1.87e+06(1.8e-08) -6.39e+05(1.3e-10)(+) -6.39e+05(1.3e-10)(+) -6.90e+05(1.5e-10)(+) sector (500-1000-900) - Justification: In the abstract, we highlighted our contributions, including algorithm development, theoretical analysis, and empirical study.Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper.• The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.• The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.• It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.Justification: Our research aligns with the ethical guidelines outlined by NeurIPS.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics.• The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction).Taking n = 4 for example, we have: Υ = {{(1,2),(3,4)}, {(1, 3), (2, 4)}, {(1, 4), (2, 3)}}. 


Decision: This paper had borderline scores, slightly positive (avg 5.25) but below the standard acceptance threshold (~5.7). That means the Area Chair took a very close look at the paper and comments.

The reviewers are mildly positive and the biggest strength seems to be that this is a pretty reasonable method to solve optimization under J-Orthogonality constraints (for reference, this means $X$ must satisfy $X^T J X$ where $J$ is block diagonal with one identity block and one negative identity block, hence "hyperbolic" style). There are not a lot of methods to solve this, and of applicable methods, the numerical results show that the proposed approach does very well.

That said, the reviewers raise two main issues:
1. the key algorithmic innovation seems to be mild, since much of this is similar to the work in [51]. The authors argue against this in their rebuttal, mentioning the addition of several things including an optimality condition they derived and variance reduction. They also claim that J-orthogonality constraints, not representing a compact set, are fundamentally different.  Reviewer swMz thought that this compactness argument didn't have a direct role. Ultimately, the rebuttal slightly increased one or two reviewers' scores but didn't make any reviewer extremely positive.
2. there were some issues about the implementations of the methods compared against, raised by a few reviewers. The authors have altered other methods a bit (which seems to have been in good faith, in order to make those methods perform well) but some reviewers were worried this made those methods potentially worse. As this was raised by more than one reviewer, I think it is significant. I'd suggest that the authors compare to two versions of existing methods: one implemented exactly as in the existing papers, and one tweaked.

A final point raised later in discussion by a reviewer is that there are techniques to convert many constrained problems into unconstrained problems, and that this should be discussed (in particular, cf. Xiao, Nachuan, Xin Liu, and Kim-Chuan Toh. "Dissolving constraints for Riemannian optimization." Mathematics of Operations Research 49.1 (2024): 366-397. )

Adding my own thoughts, I think that these hyperbolic style constraints are interesting, and I appreciate the applications listed at the top of page 2 (e.g., to the hyperbolic eigenvalue problem, ultrahyperbolic knowledge graph embedding, ...), but ultimately this is a rather specialized topic. That is fine, but I think with such a specialized topic, the bar for publication should be just a little bit higher.

Thus in the end, I am unfortunately not recommending this paper for NeurIPS this year, given that the topic is so specialized and the borderline ratings from the reviewers.

Reviews: review: paper_summary: paper_summary: This paper proposes a block coordinate descent method for solving optimization problems with J-orthogonality constraints. Several variants of the method are introduced within this framework, and convergence results are established. Extensive numerical results are also presented to demonstrate the efficiency of the proposed methods. However, I have some concerns regarding the novelty of this paper as well as the numerical results.
strength_weakness: strength_weakness: 1. My major concern is that the novelty of this paper might be insufficient since the row-based approach is very similar to that in [51], even though the two papers tackle different problems.

2.  For the numerical results shown in Table 1, the proposed method fails to return a feasible solution for some instances, such as randn(10-10-5) and w1a (2470-290-145), as well as some other instances in the appendix. This is strange since the paper describes a BCD-type method, which should always return a feasible solution.

3. The information of the reference [31] might be incorrect. 

4. For the GS-JOBCD method, there are two options for choosing $Q$, whereas J-JOBCD only has one option. The authors should provide an explanation for this difference.

5. The presentation could be further improved. Here are a few examples: the formulation of $P_i$  after equation (12) could be simplified by removing the notation $\mathrm{mat}$; it is unclear if the requirement on $\underline{Q}$ in equation (4) is sufficient to guarantee convergence (probably not, since $\underline{Q} = 0$  also satisfies this condition).
questions: questions: 1. How to perform and choose the parameters in ADMM for solving this problem? It seems that the corresponding method lacks a convergence guarantee.

2. The authors should provide more details on the implementation of J-JOBCD and GS-JOBCD. Did you use parallel techniques in J-JOBCD? What would happen if these techniques were not used? Should you compare the iteration numbers of the methods? How do you choose the parameter in equation (12)? Moreover, if J-JOBCD consistently performs much better than GS-JOBCD, it seems unnecessary to introduce GS-JOBCD.
limitations: limitations: At the beginning of the paper, the authors claim that equation (2) can imply 
$\\|\nabla f_i(X) - \nabla_i f(X^+)\\| \leq L_f \\|X - X^+\\|$, which is incorrect. Note that the converse is correct. The other assumptions in Assumptions 4.1 and 4.2 essentially assume the compactness of the iterates.
paper_summary: paper_summary: The paper proposes two block coordinate descent methods for minimization of a finite-sum subject to the J-Orthogonality constraints — one based on Gauss-Seidel strategy, the other based on variance reduction and Jacobi strategy. The convergence is proved, with a global convergence rate of O(N/\epsilon) and O(\sqrt{N}/\epsilon) respectively, and a local convergence rate that depends on the desingularization in the KL-condition assumption.
strength_weakness: strength_weakness: The paper is relatively dense, and I find it a bit hard to keep track of all the terms introduced. For instance, the parameter theta is used in the algorithms but I’m not sure where it is introduced; in Assumption 4.8 KL function is mentioned but it’s not defined…
questions: questions: The blocks in the proposed JOBCD algorithms consist of just 2 indices. I’m wondering if there is benefit in updating more than 2 indices in each iteration? 

I might miss something while reading the paper… what is the parameter theta in the algorithm 1?
limitations: limitations: Yes, the paper discusses the assumptions of the theorems.
paper_summary: paper_summary: This paper proposes two Block Coordinate gradient descent methods(BCD) for solving J-orthogonal constrained problem. One is Gauss-Seidel type, the other one is Jocobi type as well as addressing finite sum problem using variance reduction strategies. Convergence guarantees are proved with KL conditions. Numerical experiments show the advantages of the  proposed algorithm.
strength_weakness: strength_weakness: 1. This paper is based on the paper " [51] Ganzhao Yuan. A block coordinate descent method for nonsmooth composite optimization
under orthogonality constraints. ArXiv, abs/2304.03641, 2023."
The main difference is the constraint in this paper becomes J-orthogonality constraint. However, the framework follows almost the same as [51]. The authors should highlight the novelty of the algorithm or difficulty in the extension.

2. The authors of the reference [31] may be wrong. Besides, the UMCM algorithm in [31] solves orthogonal constrained problem. Is there any difference in implmenting in solving J-orthogonality problem? The objective value of UMCM is far from the JOBCD method. I'm curious about the reasons.

3. In numerical experiment, how do you select subset from the dataset, see line 309.
questions: questions: 1. Why cannot we design variance reduced algorithm for Gauss-Seidel type BCD?
2. line 323: it claims that the lower objective achieved by GS-JOBCD can be explained by the stronger stationary point result. However, it is based on KL condition. I think it is not the main reason?
3. line 324: infinite sum?
limitations: limitations: N/A
paper_summary: paper_summary: The paper introduces the JOBCD (J-Orthogonal Block Coordinate Descent) algorithm, a novel method designed to tackle optimization problems under J-orthogonality constraints. JOBCD includes two variants: GS-JOBCD (Gauss-Seidel strategy) and VR-J-JOBCD (Jacobi strategy with variance reduction). Theoretical analyses establish the algorithms' complexity and convergence, while extensive experiments show JOBCD's superior performance compared to state-of-the-art methods in various applications.
strength_weakness: strength_weakness: Some proofs for Section 4 are hard to follow.
questions: questions: 1. Do the optimization problems in (S3) of Algorithm 1 and Algorithm 2 require an especially exact solution? What happens if they are solved approximately?
2. Line 34: The notation $\mathcal{Q}(X^+; X)$ is undefined.
3. Line 109: Should the notation $\mathrm{B}$ be replaced by $\mathrm{B^t}$?
4. LIne 124, formula (4): Should $I_2$ be replaced by $I_4$?
limitations: limitations: NA

Comments: title: responce
Implementation of UMCM algorithm for Optimization Problem under J-Orthogonality Constraints
Implementation of ADMM Algorithm for Optimizaton Problem under J-Orthogonality Constraints
comment: **Question 11. Since the orthogonality constraint is compact, we can directly obtain the closed-form solution for the orthogonal tangent space." I don't think so. The closed-form solution is not due to the compactness of the orthogonal constraint.**

Response: It should be "For the orthogonality constraint which is compact, we can directly obtain the closed-form solution for the orthogonal tangent space". 

**Question 12. By the way, which Riemannian metric did you use to define the Riemannian gradient? If the metric is induced from the Euclidean one, then it is the projection of the Euclidean gradient onto the tangent space, for which you have the expression of the Riemannian gradient in Lemma E.3.**

Response:

- Projecting the Euclidean gradient onto the tangent space requires solving the following optimization problem: $\min_Y \|\|\mathbf{Y} - \mathbf{G}\|\|_F^2, \operatorname{s.t.} \mathbf{X}^\top \mathbf{J}\mathbf{Y} + \mathbf{Y}^\top \mathbf{J}\mathbf{X} = 0$. Although this problem involves linear constraints, obtaining a closed-form solution is challenging, so we chose to avoid solving it. In other words, we did not use the projection metric.

- Following the strategy outlined by [Wen2013], we derived the first-order optimality condition $\mathcal{J}(\mathbf{X}) = 0$ for any feasible solution $\mathbf{X}$ to the optimization problem, using the first-order information of the objective function and the symmetric property of the multiplier. We refer to $\mathcal{J}(\mathbf{X})$ as the Riemannian gradient.

- It should be noted that our JOBCD algorithm does not use the Riemannian gradient $\mathcal{J}(\mathbf{X})$, as we rely solely on the Euclidean gradient. The Riemannian gradient is employed only for theoretical analysis.


[Wen2013] Wen, Z., & Yin, W. (2013). A feasible method for optimization with orthogonality constraints. Mathematical Programming, 142(1), 397-434.
**Question 10. Regarding the UMCM method, the original algorithm should be implemented exactly, including the quadratic term, regardless of whether the subproblem is well defined or not. Deviating from this makes it a different algorithm altogether. I checked the code and noticed that the step size decreases by half every two iterations, which would degrade its performance. Since UMCM is not a stochastic algorithm, a diminishing step size is not necessary.**

Response:

- In the supplementary experiments, we essentially used a fixed step size of $10^{-5}$ for **both methods**, with and without the quadratic term, which is a relatively small step size.

- During the experiment design, we tested various step sizes, including $\lbrace 10^{-3}, 10^{-5}, 10^{-7} \rbrace$, and found that $10^{-5}$ generally yielded better results.

- One possible reason for requiring such a small step size is that the penalty term in the objective function is a **quartic function**, which is highly sensitive to the step size scale, making a very small step size necessary.

- In the actual experiments, we adopted a dynamic step size, gradually decreasing to $10^{-5}$ to achieve better outcomes.
"Since the orthogonal constraint is compact, we can directly obtain the closed-form solution for the orthogonal tangent space."
I don't think so. The closed-form solution is not due to the compactness of the orthogonal constraint. The tangent space has a nice decomposition formula; see [1]. Therefore, we can get the projection solution exactly. Indeed, as the tangent space is a linear subspace, as stated in your Lemma A.3, one can solve the projection problem using the standard projection technique onto a linear constraint.

By the way, which Riemannian metric did you use to define the Riemannian gradient? If the metric is induced from the Euclidean one, then it is the projection of the Euclidean gradient onto the tangent space, for which you have the expression of the Riemannian gradient in Lemma E.3.
Regarding the UMCM method, the original algorithm should be implemented exactly, including the quadratic term, regardless of whether the subproblem is well defined or not. Deviating from this makes it a different algorithm altogether. I checked the code and noticed that the step size decreases by half every two iterations, which would degrade its performance. Since UMCM is not a stochastic algorithm, a diminishing step size is not necessary.
**Question 9.  The main difference between J-orthogonality and orthogonality lies in boundedness. Could you please indicate where the boundedness makes your analysis unique compared to [51]?**

Response: To handle the J-orthogonal constraint, we introduce additional novel strategies in the proof of convergence. We illustrate these with examples:

1. **Lemma E.5:** Based on the KL assumption, the key proof step involves deriving the inequality relationship between $\operatorname{dist}(0, \nabla f^\circ(\mathbf{X}))$ and $\operatorname{dist}(0, \nabla_{\mathcal{J} }f(\mathbf{X}))$ by solving the optimal solution in Lemma E.4 regarding the tangent space projection:$
   \bar{\mathbf{Y}} = \underset{\mathbf{Y} \in T_{\mathbf{X}} \mathcal{J}}{\operatorname{argmin}} \, h(\mathbf{Y}) = \frac{1}{2} \|\|\mathbf{Y} - \mathbf{G}\|\|_F^2$. Since the orthogonal constraint is compact, we can directly obtain the closed-form solution for the orthogonal tangent space. However, for the J-orthogonal matrix, we cannot obtain the optimal solution. Therefore, in our proof, we establish an inequality based on the necessary (not sufficient) condition for the optimal solution: 
   $
   h(\bar{\mathbf{Y}}) \leq h(\mathbf{G} - \mathbf{J} \mathbf{X} \mathbf{G}^\top \mathbf{X} \mathbf{J})
   $.
 
2. **Lemma E.2:** In proving the "Riemannian Gradient Lower Bound for the Iterates Gap", we need to rely on the first-order optimality condition of the optimization problem under J-orthogonality constraints, constructed for the first time in Lemma 3.1:
   $
   0 = \nabla_{\mathcal{J}} f(\bar{\mathbf{X}}) = \nabla f(\bar{\mathbf{X}}) - \mathbf{J} \bar{\mathbf{X}} [\nabla f(\bar{\mathbf{X}})]^\top \bar{\mathbf{X}} \mathbf{J},
   $
   to complete the proof. This results in a completely different key coefficient, $\phi$, for the lemma.

3. Similarly, in proving **Lemma E.3** regarding $\operatorname{dist}(0, \nabla_{\mathcal{J}} f(\mathbf{X}))$ and $\nabla_{\mathcal{J}} \mathcal{T}(\mathbf{I}_2;\mathbf{X}^t, B)$, there are also many novel aspects, with significant differences in the proof process and the key coefficients in the results.

Moreover, the unbounded nature of optimization problem under J-orthogonality constraints adds significant complexity to the proofs related to variance reduction stochastic strategies and parallel strategies.
**Question 8.  In [31], the penalty function used is the augmented Lagrangian, whereas in this response, the Lagrangian is used.**

Response: 

-  The unconstrained objective function obtained using the multiplier correction strategy is equivalent and well-defined whether or not it includes quadratic terms. We have chosen the simpler version.

- Since we have already compared with the ADMM method based on the augmented Lagrangian function, using this more straightforward approach also helps to distinguish it from ADMM. 

- Now, we illustrate that the equivalent minimization function (without the quadratic term) is also well-defined for the HEVP problem. We have:

$\min \limits_{\mathbf{X} \in \mathbb{R}^{n \times n}} f(\mathbf{X}) - \tfrac{1}{2} \langle \mathbf{J}\mathbf{X}^\top \nabla f(\mathbf{X}), \mathbf{X}^\top\mathbf{J}\mathbf{X} - \mathbf{J} \rangle,$

where $ f(\mathbf{X}) = \text{tr}(\mathbf{X}^\top\mathbf{C}\mathbf{X}) $ and $\mathbf{C} = -\mathbf{D}^\top\mathbf{D}$.

In the above expression, the function $ f(\mathbf{X}) $ is quadratic, and the Lagrangian term contains a quartic function $-\langle \mathbf{J}\mathbf{X}^\top\mathbf{C}\mathbf{X}, \mathbf{X}^\top\mathbf{J}\mathbf{X} \rangle$. Now, we show that this expression is always non-negative:

$-\langle \mathbf{J}\mathbf{X}^\top\mathbf{C}\mathbf{X}, \mathbf{X}^\top\mathbf{J}\mathbf{X} \rangle = -\langle \mathbf{J}\mathbf{X}^\top\mathbf{C}\mathbf{X}, \mathbf{J} \rangle = -\langle \mathbf{X}^\top\mathbf{C}\mathbf{X}, \mathbf{I} \rangle = -\text{tr}(\mathbf{X}^\top\mathbf{C}\mathbf{X}).$

Since $\mathbf{C} = -\mathbf{D}^\top\mathbf{D}$, the term $-\text{tr}(\mathbf{X}^\top\mathbf{C}\mathbf{X})$ is always non-negative. Therefore, we conclude that this quartic term is always non-negative and dominates the other lower-order terms, ensuring that the objective function is well-defined.

- Upon request, we also report the results for the UMCM algorithm including the quadratic term (i.e., using the augmented Lagrangian function). Notably, we consider the following minimization problem:

$\min \limits_{\mathbf{X} \in \mathbb{R}^{n \times n}} f(\mathbf{X}) - \tfrac{1}{2} \langle \mathbf{J}\mathbf{X}^\top \nabla f(\mathbf{X}), \mathbf{X}^\top\mathbf{J}\mathbf{X} - \mathbf{J} \rangle + \tfrac{\beta}{2} \|\| \mathbf{X}^\top\mathbf{J}\mathbf{X} - \mathbf{J} \|\|^2_\mathsf{F}.$

Please refer to the table below. These results will be included in the revised manuscript.


| HEVP-datasetname | cifar          | CnnCaltech     | gisette        | mnist          | randn10        |
|------------------|----------------|----------------|----------------|----------------|----------------|
| (m-n-p)          | (1000-100-50)  | (2000-1000-500) | (3000-1000-500) | (1000-780-390) | (10-10-5)      |
| Lagrangian       | -4.81e+02(1.0e-06) | -8.61e+01(8.6e-07) | -1.70e+06(8.1e-07) | -2.79e+04(1.0e-06) | -3.38e-01(1.0e-06) |
| Augmented Lagrangian | -4.86e+02(1.0e-06) | -7.33e+01(3.4e-07) | -1.33e+06(5.0e-07) | -2.80e+04(1.0e-06) | -1.41e+00(1.0e-06) |

| HEVP-datasetname | randn100       | randn1000      | sector         | TDT2           | w1a            |
|------------------|----------------|----------------|----------------|----------------|----------------|
| (m-n-p)          | (100-100-50)   | (1000-1000-500) | (500-1000-500)  | (1000-1000-500) | (2470-290-145) |
| Lagrangian       | -1.13e+02(1.0e-06) | -1.28e+04(6.8e-08) | -1.39e+03(6.8e-07) | -1.73e+06(6.8e-07) | -2.89e+02(1.0e-06) |
| Augmented Lagrangian | -1.42e+02(1.0e-06) | -1.11e+04(4.9e-08) | -1.11e+03(4.2e-07) | -1.38e+06(4.3e-07) | -3.08e+02(1.0e-06) |

*Supplementary experiments for **HEVP** (limited to 30s). The data in the cell stand for $f(\mathbf{X}^t) - f(\mathbf{X}^0)$ of **UMCM** with Lagrangian and Augmented Lagrangian ($\beta = 10$) methods. The value in parentheses represents $\tfrac{1}{n^2} \sum_{ij}^n |\mathbf{X}^{\top}\mathbf{J}\mathbf{X}-\mathbf{J}|_{ij}$.*
Thank you to the authors for the clarification. Some of my concerns have been well addressed. Therefore, I would like to increase my score to 5.
Thank you for the clarification. However, I noticed that in [31], the penalty function used is the augmented Lagrangian, whereas in this response, the Lagrangian is used. Additionally, you mentioned that the main difference between J-orthogonality and orthogonality lies in boundedness. Could you please indicate where the boundedness makes your analysis unique compared to [51]?
Thank you to the authors for their responses. I will maintain my score.
Thank the authors for the clarifications! In the revised version, please be sure to properly introduce parameters/variables/terms before using them (such as theta and KL function). I will keep my score.
We consider minimizing the following differentiable function under J-orthogonality constraints:

$\min \limits_{\mathbf{X} \in \mathbb{R}^{n \times n}}  f\mathbf{(X)} , \text { s.t. }  \mathbf{X}^{\top} \mathbf{J} \mathbf{X}=\mathbf{J}$. 

We derive the Lagrangian function of the above problem with $\Lambda \in \mathbb{R}^{n \times n}$:

$\mathcal{L}(\mathbf{X}, \Lambda) = f(\mathbf{X})-\tfrac{1}{2} \langle \Lambda, \mathbf{X}^\top\mathbf{J}\mathbf{X} - \mathbf{J} \rangle. $

Setting the gradient of $\mathcal{L}(\mathbf{X}, \Lambda)$ w.r.t. $\mathbf{X}$ to zero yields:

$\nabla f(\mathbf{X}) - \mathbf{J} \mathbf{X} \Lambda = 0. $

Multiplying both sides by $\mathbf{X}^\top$ and using the fact that $\mathbf{X}^\top \mathbf{J}\mathbf{X}=\mathbf{J}$, we have $\mathbf{J} \Lambda = \mathbf{X}^\top \nabla f(\mathbf{X})$. Multiplying both sides by $\mathbf{J}^\top$ and using $\mathbf{J}^\top \mathbf{J}=\mathbf{I}$, we have $\Lambda=\mathbf{J}\mathbf{X}^\top \nabla f(\mathbf{X})$. Thus, we obtain the following equivalent unconstrained optimization problem: 

$\min \limits_{\mathbf{X} \in \mathbb{R}^{n \times n}} f\mathbf{(X)}-\tfrac{1}{2} \langle \mathbf{J}\mathbf{X}^\top \nabla f(\mathbf{X}), \mathbf{X}^\top\mathbf{J}\mathbf{X} - \mathbf{J} \rangle. $

Finally, we solve it using a gradient-based approach, specifically employing the Adagrad optimizer built into PyTorch in our paper.
We consider minimizing the following differentiable function subject to J-orthogonality constraints:

$\min \limits_{ X \in \mathbb{R}^{n \times n}} f(X),
\text { s.t.}  X^{\top}  J   X= J $.

Defining $ Y = J   X \in \mathbb{R}^{n \times n}$, we have:
$\min \limits_{ X , Y  \in \mathbb{R}^{n \times n}}  f(X),
\text { s.t. }  X ^{\top}  Y = J ,
  Y = J   X $.

Introducing Lagrangian multipliers $ Z  \in \mathbb{R}^{n \times n}$ and $ W  \in \mathbb{R}^{n \times n}$, and the penality variable $\beta \in \mathbb{R}$, we have the following argumented Lagrangian function:

$\mathcal{L}( X , Y ; Z , W ) = f( X )+\langle X ^\top Y - J , Z \rangle+\langle J  X - Y , W \rangle+\frac{\beta}{2}\|\| X ^\top Y - J \|\|_F^2+\frac{\beta}{2}\|\| J  X - Y \|\|_F^2$

Suppose $f(X)$ is $l$-Lipschitz gradient continuous: $f(X) \leqslant f( X ^t)+ \langle  X - X ^t, \nabla f ( X ^t ) \rangle+\frac{l}{2}\| \| X - X ^t \|\|_{F}^2$.

We consider minimizing $\mathcal{L}(X,Y ;Z,W )$ with respect to $X$. We derive the following majorization function $Q^t(X)$ of $\mathcal{L}(X,Y^t;Z^t,W^t)$ at $X^t$:

$\mathcal{L}(X, Y^t ; Z^t, W^t) \leq  Q^t(X) \triangleq f (\mathbf{X}^t)+ \langle  X - X ^t, \nabla f (X^t ) \rangle+\frac{l}{2} \|\| X - X ^t \|\|_{F}^2+ \langle X^\top Y^t - J , Z^t \rangle+\langle J  X - Y^t, W^t \rangle+\frac{\beta}{2}\|\| X^\top Y^t - J \|\|_F^2+\frac{\beta}{2}\| \|J  X - Y^t \|\|_F^2 $

We solve the following subproblem to update $ X ^{t+1}$ and $ Y ^{t+1}$ alternately:

$ X ^{t+1} = \operatorname{arg} \operatorname{min}_ X  Q^t(X)$


$ Y ^{t+1} = \operatorname{arg} \operatorname{min}_ Y P^t(Y) \triangleq \mathcal{L}( X^t , Y ; Z^t , W^t ) $


$ Z ^{t+1} =  Z ^t +\beta\cdot({ X ^{t+1}}^\top Y ^{t+1}- J ) $

$ W ^{t+1} =  W ^t +\beta\cdot( J { X ^{t+1}}- Y ^{t+1}) $

Using the first-order optimality conditions of $\operatorname{min}_ X  Q^t(X)$ and $\operatorname{min}_ Y P^t(Y)$, we get the updated formula for $ X ^{t+1}$ and $ Y ^{t+1}$:

$ X ^{t+1} = -(l\mathbf{I}+\beta( Y ^t{ Y ^t}^\top+\mathbf{I}))^{-1}(\nabla f({ X ^t}^\top)-l X ^t+ Y ^t{ X ^t}^\top+ J  W ^t-\beta Y ^t J -\beta J  Y ^t) $

$ Y ^{t+1} = -(\beta( X ^{t+1}{ X ^{t+1}}^\top+\mathbf{I}))^{-1}( X ^{t+1} Z ^t- J  W ^t-\beta X ^{t+1} J -\beta J  X ^{t+1}) $

In the experiment, we used the following strategy to adjust $\lambda$: If $t$\%$2 ==0$ and $\lambda \leq 1e9$, then $\lambda$ is doubled.
**Question 8. Did you use parallel techniques in J-JOBCD?**

Response: Yes, we utilized parallel techniques in J-JOBCD by expanding the matrix dimensions and defining the variables as tensors in  PyTorch to facilitate parallel computations.


**Question 9. What would happen if these techniques were not used?**

Response: 

1. Without using parallel techniques, the method would degrade into a sequential BCD method. The parallelization strategy in this paper is an effort to accelerate the BCD framework by fully utilizing modern parallel architectures.

2. We could consider using other parallelization techniques for experiments, such as more efficient multi-core parallel computing or efficient GPU-based parallel computing. It is worth noting that using PyTorch to define variables as tensors for parallel computations is very efficient for small-dimensional computations, as it avoids the communication costs between different machines.

**Question 10. Should you compare the iteration numbers of the methods?**

Response: Thank you for your suggestion. **We have conducted relevant experiments, as shown in Table 3 and Figure 2 of the attached PDF file**. The JOBCD algorithm consistently achieves better optimization results in fewer iterations.


**Question 11. How do you choose the parameter in equation (12)?**

Response: In the experiments of the paper, we set $\theta = 10^{-6}$ and let $\zeta$ be the Lipschitz constant of the blocks of different functions. Specifically, in the HEVP experiments, we used the maximum absolute value of the chosen blocks in the Hessian matrix $C$ as the block Lipschitz constant. In the HSPP experiments, we heuristically specified the Lipschitz constants empirically for different datasets and selected the most effective constant. We chose within the range of $\lbrace 100, 500, 1000 \rbrace$, and generally, 500 yielded the best results.

**Question 12. Moreover, if J-JOBCD consistently performs much better than GS-JOBCD, it seems unnecessary to introduce GS-JOBCD.**

Response: We could avoid GS-JOBCD and directly present the J-JOBCD algorithm. However, doing so would not illustrate the practical impact of parallelization and variance reduction techniques on the algorithm's convergence. 

To better demonstrate the influence of these modern techniques, we described GS-JOBCD and J-JOBCD separately in the writing. 

Additionally, GS-JOBCD can serve as a basic deterministic algorithmic framework for further designs.

**Question 12. At the beginning of the paper, the authors claim that equation (2) can imply $|| \nabla f_i(X)  -\nabla f_i(X^+) ||_F\leq L_f ||X-X^+||_F$, which is incorrect. Note that the converse is correct.**

Response: We will explicitly assume $|| \nabla f_i(X)  -\nabla f_i(X^+) ||_F\leq L_f ||X-X^+||_F$ for some constant $L_f$. Thank you for pointing it out.