OPTIMIZATION OVER THE INTERSECTION OF MANIFOLDS∗
arXiv:2605.22736v1 [math.OC] 21 May 2026
YAN YANG† , BIN GAO‡ , AND YA-XIANG YUAN‡ Abstract. Optimization over the intersection of two manifolds arises in a broad range of applications, but is hindered by the coupled geometry of the feasible region. In this paper, we prove that the regularities—clean intersection and intrinsic transversality—are equivalent, which yields a tractable projection onto the tangent space of the intersection. Therefore, we propose a geometric method that employs a retraction on only one manifold and updates the iterate along two orthogonal directions. Specifically, the iterates stay on one manifold, and the two directions are responsible for asymptotically approaching the other manifold and decreasing the objective function, respectively. Under intrinsic transversality, we derive the convergence rate for both the feasibility and optimality measures, and show that every accumulation point is first-order stationary. Numerical experiments on problems stemming from sparse and low-rank optimization, including fitting spherical data, approximating hyperbolic embeddings on real data, and computing compressed modes, demonstrate the effectiveness of the proposed method. Key words. Manifold intersection, tangent space, orthogonal directions, intrinsic transversality, clean intersection, Riemannian optimization MSC codes. 65K05, 90C30, 90C46
1. Introduction. In this work, we consider the following optimization problem over the intersection of two constraint sets: min
f (X)
s. t.
h(X) = 0, X ∈ M,
X∈E
(P)
where E denotes a finite-dimensional Euclidean space, accommodating Rm , Rm×n , or Rn1 ×···×nd , and M ⊆ E is a smooth submanifold. The objective function f : E → R and the constraint-defining map h : E → Rq are both smooth. We denote the zero level set of h by H := {X ∈ E | h(X) = 0} , and thus the feasible region is the intersection H ∩ M. Throughout this paper, we impose the following assumption. Assumption 1. There exists an open neighborhood K ⊆ E of H such that the differential DhX : E → Rq has full rank q for all X ∈ K.
Assumption 1 implies that H is a smooth manifold in E; see [32, Corollary 5.14]. In the vanilla scenario h(·) ≡ 0, i.e., H = E, (P) reduces to an unconstrained optimization problem on the smooth manifold M, for which a variety of algorithms— including Riemannian gradient descent and trust-region methods—are well established; see [3, 12]. However, once the additional constraint X ∈ H is non-trivial, ∗ Submitted to the editors DATE.
Funding: This work was supported by the National Key R&D Program of China (grant 2023YFA1009300). BG and YY were supported by the National Natural Science Foundation of China (grant No. 12288201). † State Key Laboratory of Mathematical Sciences, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, and University of Chinese Academy of Sciences, China ([email protected]). ‡ State Key Laboratory of Mathematical Sciences, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, China ([email protected], [email protected]). 1
2
Y. YANG, B. GAO, AND Y.-X. YUAN
dealing with the intersection H ∩ M presents several challenges. First, the intersection H ∩ M does not necessarily constitute a smooth manifold, which impedes the direct application of existing Riemannian optimization algorithms. Second, the geometry of H ∩ M is more intricate than that of M or H alone—for instance, unclear characterization of the tangent cone to the intersection hinders the construction of effective search directions. Third, projections onto H ∩ M are generally unavailable in closed form, making it difficult to preserve the feasibility of the iterates. 1.1. Motivation and applications. The formulation (P) encompasses a broad range of problems where a manifold constraint X ∈ M is coupled with additional structured requirements h(X) = 0. We outline several representative applications. Table 1 Instances of manifold intersections arising in sparse optimization, low-rank matrix and tensor optimization. The intersection geometry of the specific (M, H) pairs is summarized. Specifically, all the intersections satisfy the intrinsic transversality.
Sparse
Low-rank
General
Manifold
Level set of h
Cs
{X ∈ Rm | ∥X∥2F = 1} n×p
Intersection geometry
⊤
| X X = Ip }
[10]
Cs
{X ∈ R
Mr
{X ∈ Rm×n | A(X) = b}
[36, 64]
Mr
H is orthogonally invariant
[65, 64]
Mr
H is hyperbolic (7.2)
Appendix A
Sr (n)
{X ∈ S(n) | ∥X∥2F = 1}
[16, 37, 64]
S+ r (n)
{X ∈ S(n) | A(X) = b}
[33, 64]
n1 ×n2 ×···×nd
Mtt r
| ∥X∥2F = 1}
{X ∈ R
M
H satisfies Assumption 1
[18]
[49, 64] Intrinsic transversality
Let ∥ · ∥0 denote the cardinality of an element. The sparsity set M = Cs := {X ∈ E | ∥X∥0 = s} combined with normalization or orthogonality constraints appears in several scenarios. When E = Rm is the vector space, Beck and Hallak [10] characterized the projection onto Cs ∩ {X ∈ Rm | ∥X∥2F = 1}, with applications in the fields of genetics and finance. Moreover, when E = Rn×p , the intersection of Cs and the Stiefel manifold St(n, p) := {X ∈ Rn×p | X ⊤ X = Ip } underlies sparse principal component analysis [19] and the geometry of St(n, p) ∩ Cs has been recently studied in [18]. Another important class of instances concerns the fixed-rank manifold M = Mr := {X ∈ Rm×n | rank(X) = r}. Specifically, Cason et al. [16] derived the tangent cone to H ∩ Mr with H as the Frobenius sphere, applied to the approximation of graph similarity matrices. Li and Luo [36] subsequently obtained the normal cone to H ∩ Mr for H as an affine manifold. More generally, Yang et al. [65] characterized the tangent and normal cones to H ∩ Mr when h is orthogonally invariant, i.e., h(X) = h(XQ) for all orthogonal Q. In addition, the hyperbolic constraint, arising from low-rank compression of hyperbolic embeddings for hierarchical data [29], was treated in [64]. In the symmetric setting, Li et al. [37] studied the geometry of Sr (n) := {X ∈ S(n) | rank(X) = r} intersected with spectral constraints, where S(n) := {X ∈ Rn×n | X ⊤ = X}. Levin et al. [33] analyzed the positive semidefinite counterparts, denoted by H ∩ S+ r (n) for some H. The results were extended to ten-
OPTIMIZATION OVER THE INTERSECTION OF MANIFOLDS
3
sors: the set Mtt r of low-rank tensors in tensor-train format coupled with a sphere was investigated in [49] with applications to quantum physics. A unified analysis on the geometry of H ∩ M covering the above low-rank instances was developed in [64]. We refer the reader to Table 1 for a summary. 1.2. Related work. We begin with M = E and the role of decomposing search directions into orthogonal components for tackling the equality constraint h(X) = 0. Direction decomposition for equality-constrained problems. If the region M is the whole Euclidean space E, (P) reduces to an optimization problem solely subject to the equality constraint h(X) = 0. A representative principle is to decompose the update direction into tangent and normal components, responsible for decreasing the objective f and for encouraging the feasibility X ∈ H, respectively. This idea was first formalized by [52, 53], and then Frost [24] proposed the corrective gradient projection method, which realized a normal direction as the correction to the tangent one. Additionally, second-order information was exploited in the so-called null-space methods [43, 66], which enhanced the convergence results. More relevant to our work, the landing algorithm was proposed for optimization problems with orthogonality constraints [1, 25, 2], getting rid of computationally expensive retractions adopted by Riemannian optimization methods [3, 12]. Then, Schechtman et al. [54] extended the technique to general equality constraints, accommodating stochastic oracles, and Vary et al. [59] adapted the algorithm for problems over the random generalized Stiefel manifold, which was further extended to distributed optimization [57]. Subsequent work [56, 26] incorporated backtracking line search into the landing method. More recently, Xiong et al. [62] designed a second-order landing algorithm to achieve locally quadratic convergence. Goyens and Feppon [26] unveiled the relationship between the landing algorithm and several classical optimization methods including the sequential quadratic programming method and the augmented Lagrangian method [44]. When the manifold constraint X ∈ M is non-trivial, i.e., M ̸= E, finding feasible points in H ∩ M is itself a challenging problem, and thus we review the theory of the method of alternating projections (MAP) [41]. Intersection condition and alternating projection. Given two general closed sets H and M in a Euclidean space, the feasibility problem seeks a point X ∗ ∈ H ∩ M. The method of alternating projections generates iterates by Xk+1 ∈ PM (PH (Xk )), where P denotes the projector onto a closed set. When H and M are nonconvex, establishing local linear convergence requires appropriate regularity conditions of the intersection H ∩ M, often realized as a separation property of the limiting normal cones. The conditions proposed in existing work are summarized in Figure 1. Specifically, Lewis and Malick [35] first established local linear convergence of the MAP under transversality when both H and M are manifolds, which was generalized to clean intersection by Andersson and Carlsson [4]. In parallel, Lewis et al. [34] introduced linear regularity for general closed sets, and Bauschke et al. [9, 8] weakened it to the restricted regularity. The regularity conditions were further weakened to intrinsic transversality [21] and the separable condition [45], respectively. In addition, the local convergence can be preserved when the exact projections are replaced with appropriate inexact ones [22, 15, 61]. The recent work [17] reveals that the alternating projections can further induce retractions over manifold intersections. To the best of our knowledge, intrinsic transversality appears to be one of the most general conditions to guarantee local linear convergence of the MAP. To address the optimization problem (P), existing methods harness the specific structure of the constraint sets; we summarize them as follows.
4
Y. YANG, B. GAO, AND Y.-X. YUAN Transversality [35]
Clean intersection [4]
Intrinsic transversality [21]
Linear regularity [34]
Restricted regularity [9]
Separable condition [45]
Fig. 1. Development of intersection regularity conditions for local linear convergence of the MAP. An arrow from A to B indicates that A implies B.
Optimization over the intersection of sets. For a closed and convex M ⊆ E, a variety of methods have been developed for (P); see [50, 5, 60]. When a smooth manifold M is considered (generally nonconvex), a common perspective in recent literature is to cast (P) as a Riemannian optimization problem on M with nonlinear constraints h(X) = 0. The constraint qualifications and optimality conditions have been investigated in [63, 11, 6, 7]. Several Riemannian augmented Lagrangian methods were proposed [38, 67, 30, 6, 7], which handle the constraints h(X) = 0 by addressing a sequence of surrogate subproblems in the form of minX∈M f (X) + ⟨λ, h(X)⟩ + µ2 ∥h(X)∥2 . Moreover, second-order methods have also been developed for (P), including Riemannian sequential quadratic optimization methods [55, 46] and Riemannian interior point methods [31]. However, all the aforementioned Riemannian methods require solving a subproblem at each iteration. 1.3. Contributions. In this work, we aim to develop a geometric method for (P) by exploiting the structure of the manifold intersection. Under Assumption 1, the differential of h has full rank in the neighborhood K, and thus the set HX := {X̃ ∈ E | h(X̃) = h(X)}, serving as a perturbation of H, is a smooth manifold for all X ∈ K. In addition, we assume the intrinsic transversality condition (see Definition 2.2) as follows. Assumption 2. For all X ∈ K ∩ M, the manifolds HX and M are intrinsically transversal at X ∈ HX ∩ M.
To alleviate the difficulty arising from the intricate coupling of H and M in the constraints of (P), we preserve X ∈ M along the iterates via the retraction on M, in the spirit of the Riemannian optimization framework. More importantly, we interpret (P) as two sub-tasks: identifying feasible points in H∩M and decreasing the objective f . This perspective, together with the decomposition principle introduced in section 1.2, inspires us to seek two orthogonal directions in the tangent space of M that handle the two tasks respectively: a feasibility direction Gh that drives the iterates towards H ∩ M, and an optimality direction Gf that accounts for the descent of f . The resulting update rule takes the following form, Xk+1 = RM Xk (αk Gh (Xk ) + βk Gf (Xk )), where αk and βk are the step sizes. We then concentrate on constructing the two orthogonal directions tangent to M, which resorts to the intersection geometry. Central to the development are two new equivalent characterizations of intrinsic transversality. For the feasibility direction, we project the Gauss–Newton direction—an approximation of (PH (X) − X) that pushes X toward H—onto the tangent space of M. We prove in Theorem 4.4 that intrinsic transversality is equivalent to a projection-based
OPTIMIZATION OVER THE INTERSECTION OF MANIFOLDS
5
transversality condition; this equivalence ensures that the projected Gauss–Newton direction retains a sufficient tangential component, thereby providing an effective improvement on the feasibility. For the optimality direction, we project the negative gradient −∇f (X) onto the tangent cone THX ∩M (X). To compute this projection, we establish in Theorem 5.1 that, given two manifolds generally, intrinsic transversality is equivalent to clean intersection, answering an open question posed in [28, §8]. Therefore, HX ∩ M is a manifold with the following intersection rule, THX ∩M (X) = THX (X) ∩ TM (X),
(1.1)
Consequently, the above identity reveals that THX ∩M (X) is a linear subspace, and the projection onto it is characterized explicitly in Proposition 5.2. Combining the feasibility and optimality directions in the tangent space, we propose the Geometric method via Orthogonal Tangent Directions (GOTD), with the iterates staying on M, asymptotically approaching H, and decreasing f . The implementation is summarized in Algorithm 1 and is illustrated in Figure 2. Moreover, we √ prove an O(1/ K) convergence rate for both the feasibility and the optimality measures in Theorem 6.8, and show that every accumulation point is first-order stationary for (P) in Corollary 6.10, under some constraint qualifications. Numerical experiments on fitting spherical data, approximating hyperbolic embeddings on real data, and computing compressed modes demonstrate the effectiveness and efficiency of GOTD, attributable to the exploitation of the intersection geometry. Gf TM (Xk )
Xk
Gh
H
Xk+1
H∩M X⋆ M
Fig. 2. Illustration of GOTD. At Xk ∈ M, the update direction combines two orthogonal components in the tangent space TM (Xk ): the feasibility direction Gh that drives Xk towards H and the optimality direction Gf that decreases f . The Xk+1 is obtained by a retraction on M.
1.4. Organization. Section 2 presents the notation and preliminaries. We propose the main framework in section 3, and construct in sections 4–5 the feasibility and optimality directions via two equivalent characterizations of intrinsic transversality. Section 6 provides the convergence analysis. Section 7 validates the method on sparse and low-rank optimization problems, and the conclusion is drawn in section 8. 2. Notation and preliminaries. We outline the notation adopted in this paper, and then review some background in variational analysis and Riemannian optimization; see [3, 35, 51, 32, 12] for more details and references. 2.1. Notation. The diagonal matrix with entries x is written as Diag(x), and diag(X) extracts the diagonal of X as a vector. On a Euclidean space, we adopt
6
Y. YANG, B. GAO, AND Y.-X. YUAN
the Frobenius inner product ⟨X1 , X2 ⟩ := tr(X1⊤ X2 ), with the induced norm ∥X∥ := p ⟨X, X⟩. For a set X ⊆ E, the distance from Y to X is dist(Y, X ) := inf X∈X ∥X − Y ∥, and PX stands for the projection onto X . When X is additionally a smooth manifold, TX (X) refers to its tangent space at X, and any smooth map F : X1 → X2 between manifolds admits the differential DFX : TX1 (X) → TX2 (F (X)) at X. Given a matrix X of rank s, we write its singular value decomposition by X = U ΣV ⊤ with U ∈ St(m, s), Σ ∈ Rs×s , and V ∈ St(n, s); then the Moore–Penrose inverse is X † = V Σ −1 U ⊤ . Given a map F : E → E ′ between two Euclidean spaces, we use im(F ) ⊆ E ′ to denote the image and F ∗ : E ′ → E to denote the adjoint operator. The operator ⊙ denotes the Hadamard (entry-wise) product. 2.2. Preliminaries. Let X be a locally closed set in a finite-dimensional Euclidean space E, i.e., every point in X admits a closed neighborhood B ⊆ E such that B ∩ X is closed in E. The Bouligand tangent cone to X at a point X ∈ X is (2.1)
TX (X) := {η ∈ E | there exists ti → 0, such that dist(X + ti η, X ) = o(ti )} .
Taking the polar operation on TX (X) yields the Fréchet normal cone, ◦
NX (X) := (TX (X)) = {Y ∈ E | ⟨Y, η⟩ ≤ 0, for all η ∈ TX (X)} . When X is a smooth manifold, the tangent and normal cones coincide with the tangent and normal spaces, respectively. We then consider the intersection of two sets X and Y. Given X ∈ X ∩ Y, it generally holds that (2.2)
TX ∩Y (X) ⊆ TX (X) ∩ TY (X)
and
NX ∩Y (X) ⊇ NX (X) + NY (X).
If X and Y are smooth manifolds and they intersect transversally at X ∈ X ∩ Y, i.e., (2.3)
TX (X) + TY (X) = E, or equivalently, NX (X) ∩ NY (X) = {0},
then X ∩ Y is also a smooth manifold around X with (2.4)
TX ∩Y (X) = TX (X) ∩ TY (X)
and
NX ∩Y (X) = NX (X) + NY (X).
Note that the definition of transversality (2.3) depends on the ambient space E, and thus a generalization called clean intersection [27] was considered in [4] for the convergence analysis of alternating projections. Definition 2.1 (Clean intersection). Two manifolds X , Y ⊆ E intersect cleanly at Z ∈ X ∩ Y if X ∩ Y is a smooth manifold in a neighborhood B of Z and it holds that TX ∩Y (X) = TX (X) ∩ TY (X) for all X ∈ X ∩ Y ∩ B. Moreover, Drusvyatskiy et al. [21] introduced the following notion, which characterizes the intersection via pairs of nearby points. Definition 2.2 (Intrinsic transversality). Two locally closed sets X , Y ⊆ E are intrinsically transversal at Z ∈ X ∩ Y if there exists a constant κ ∈ (0, 1] such that, for all X ∈ X \ Y and Y ∈ Y \ X in a neighborhood of Z, (2.5)
max {dist(u, NY (Y )), dist(u, −NX (X))} ≥ κ,
where u = (X − Y )/∥X − Y ∥. The κ is called the constant of intrinsic transversality.
7
OPTIMIZATION OVER THE INTERSECTION OF MANIFOLDS
Intrinsic transversality reveals that the difference direction u cannot lie close to NY (Y ) and −NX (X) simultaneously, reflecting a separation property of the two cones. When X and Y are manifolds, clean intersection implies intrinsic transversality [21]. For problems constrained on a smooth manifold X , the framework of Riemannian optimization is developed, leveraging the Riemannian geometry of manifolds; see [3, 12] for an overview. To guide the movement from the current point along a tangent vector, a geometric tool retraction is introduced.S Specifically, a smooth map RX : TX → X , defined on the tangent bundle TX := X∈X TX (X), is called a retraction on the manifold X if for all X ∈ X , ξ ∈ TX (X), the curve γ(t) := RX X (tξ) satisfies X γ(0) = X and γ ′ (0) = ξ, where RX denotes the restriction of R on T (X). X X 3. A geometric framework via orthogonal tangent directions. We now present a geometric framework for problem (P). To preserve the structure of M, we adopt the following retraction-based update rule: (3.1)
Xk+1 = RM Xk (αk Gh (Xk ) + βk Gf (Xk )),
where Gh (Xk ), Gf (Xk ) ∈ TM (Xk ). Interpreting (P) as two sub-tasks—decreasing dist(·, H) for feasibility and decreasing f (·) for optimality—we then design Gh and Gf to handle the tasks respectively, inspired by the decomposition principle introduced in section 1.2. Under Assumptions 1–2 and given an iterate X ∈ K ∩ M, we recall the level set HX = {X̃ ∈ E | h(X̃) = h(X)} and denote the tangent cone by S(X) := THX ∩M (X). The following observation motivates the construction. Lemma 3.1. Given two manifolds X and Y in E with a point X ∈ X ∩ Y. Then for all d ∈ NY (X) and η ∈ E, it holds that PTX (X) (d), PTX ∩Y (X) (η) = 0. Proof. Let η̄ = PTX ∩Y (X) (η). By (2.2), we have η̄ ∈ TX ∩Y (X) ⊆ TX (X), and thus ⟨PTX (X) (d), η̄⟩ = ⟨d, η̄⟩. Similarly, the inclusion η̄ ∈ TY (X) holds, which together with d ∈ NY (X) yields ⟨d, η̄⟩ = 0.
Applying Lemma 3.1 with (X , Y) = (M, HX ), the orthogonality ⟨Gh , Gf ⟩ = 0 is guaranteed whenever the two directions take the following form, (3.2)
Gh (X) = PTM (X) (d(X)) and Gf (X) = PS(X) (η(X)),
where d(X) ∈ NHX (X) and η(X) ∈ E. Therefore, we shift our focus toward the design of d(X) and η(X) in sections 4 and 5, respectively, which mainly rely on the geometry of the manifold intersection. The overall development is illustrated in Figure 3, and the geometric framework is realized in Algorithm 1. Feasibility
Theorem 4.4
dist(·, H) ↓
Gh (X) ∈ TM (X)
orthogonal
Problem (P) Optimality f (·) ↓
Theorem 5.1
GOTD Algorithm 1
Gf (X) ∈ TM (X)
Fig. 3. Development of the GOTD algorithm: two equivalent characterizations of intrinsic transversality guide the construction of the orthogonal tangent directions Gh and Gf .
We outline some lemmas in preparation for sections 4–5; the reader is invited to refer back as appropriate.
8
Y. YANG, B. GAO, AND Y.-X. YUAN
3.1. Auxiliary results on differential manifolds. The following lemma reveals that the difference between two nearby points is nearly tangent. Lemma 3.2. Let X be a C 2 embedded submanifold of E and let Z ∈ X . Then there exist a neighborhood B of Z and a constant CX > 0 such that for all X1 , X2 ∈ X ∩ B, dist X1 − X2 , TX (X2 ) ≤ CX ∥X1 − X2 ∥2 . Proof. By the C 2 regularity of X , there exists a neighborhood B of Z and a C 2 parametrization ψ of X ∩B with ψ(0) = Z and TX (Z) = im(Dψ(0)). Write Xi = ψ(si ) for i = 1, 2. A Taylor expansion gives X1 − X2 = Dψ(s2 )[s1 − s2 ] + O(∥s1 − s2 ∥2 ). Since Dψ(s2 )[s1 −s2 ] ∈ TX (X2 ) and ∥s1 −s2 ∥ = O(∥X1 −X2 ∥), the estimate follows.
The next lemma collects properties of the projection operator near a C 2 submanifold.
Lemma 3.3. Let X be a C 2 embedded submanifold of E and let Z ∈ X . Then there exist a neighborhood B of Z, constants δ > 0 and L > 0 such that the projection PX is single-valued and L-Lipschitz on {X̃ ∈ E | dist(X̃, X ∩ B) < δ}. Moreover, for every X ∈ X ∩ B, unit vector v ∈ NX (X), and t ∈ (0, δ), one has PX (X + tv) = X. Proof. This is a direct consequence of [23, Theorem 4.8].
We then recall that a submanifold can be written as a graph over its tangent space. Lemma 3.4 (Local graph representation). Let X be a C k (k ≥ 1) embedded submanifold of E and let Z ∈ X . Let L := TX (Z). Then there exist a neighborhood UL ⊂ L of 0 and a C k map µ : UL → L⊥ with µ(0) = 0 and Dµ(0) = 0 such that X ∩ B = {Z + u + µ(u) | u ∈ UL } for some neighborhood B of Z.
Proof. Define π : E → L by π(X) := PL (X − Z). For all v ∈ TX (Z) = L, we have D(π|X )(Z)[v] = PL (v) = v, i.e., D(π|X )(Z) = idL . By the inverse function theorem, π|X is a local C k diffeomorphism near Z with π(Z) = 0. Let ϕ : UL → X be its local inverse, such that π(ϕ(u)) = u for all u ∈ UL . Define µ(u) := ϕ(u) − Z − u. Then µ(u) ∈ L⊥ since PL (µ(u)) = π(ϕ(u)) − u = 0, and µ is C k with µ(0) = ϕ(0) − Z = 0. For Dµ(0): since ϕ maps into X and ϕ(0) = Z, the image of Dϕ(0) lies in TX (Z) = L; differentiating π ◦ ϕ = id gives PL ◦ Dϕ(0) = idL , and thus Dϕ(0) = idL . It follows that Dµ(0)[v] = Dϕ(0)[v] − v = 0 for all v ∈ L. 3.2. Regularities preserved under diffeomorphisms. We show in the following lemma that clean intersection is preserved under local diffeomorphisms. Lemma 3.5. Let X , Y ⊆ E be C 1 embedded submanifolds that intersect cleanly at Z ∈ X ∩ Y, and let Ψ : E → E ′ be a C 1 diffeomorphism defined near Z. Then Ψ(X ) and Ψ(Y) intersect cleanly at Ψ(Z) in E ′ .
Proof. Since Ψ is a diffeomorphism, Ψ(X ∩ Y) = Ψ(X ) ∩ Ψ(Y) locally, which is a C 1 submanifold. Moreover, the chain rule gives TΨ(X ) (Ψ(Z)) = DΨ(Z)[TX (Z)] and similarly for Y and X ∩ Y. Since DΨ(Z) is an isomorphism, we have DΨ(Z)[TX (Z)] ∩ DΨ(Z)[TY (Z)] = DΨ(Z)[TX (Z) ∩ TY (Z)] = DΨ(Z)[TX ∩Y (Z)].
The same argument holds for points around Z. The next lemma, noted as an exercise in [21, §3], reveals that intrinsic transversality is preserved under local diffeomorphisms; we restate it below for completeness. Lemma 3.6. Let X , Y ⊆ E be C 1 embedded submanifolds that are intrinsically transversal at Z ∈ X ∩ Y, and let Ψ : E → E ′ be a C 1 local diffeomorphism defined near Z. Then Ψ(X ) and Ψ(Y) are intrinsically transversal at Ψ(Z) in E ′ .
OPTIMIZATION OVER THE INTERSECTION OF MANIFOLDS
9
4. Projection-based transversality and feasibility direction. By (3.2), designing the feasibility direction Gh reduces to specifying d(X). Since the goal of Gh is to drive the iterate X ∈ M towards H, a reasonable candidate is along the projection residual r(X) := PH (X) − X, which points from X towards H. However, recalling that d(X) is required to lie in NHX (X) to guarantee the orthogonality of Gh and Gf , we instead approximate r(X) by the Gauss–Newton direction, which solves the linearized least-squares problem mind∈NHX (X) ∥h(X) + DhX (d)∥2 : (4.1)
d(X) = −Dh∗X DhX Dh∗X
−1
h(X).
Note that the direction d(X) lies in NHX (X) = im(Dh∗X ), and it approximates the projection residual r(X) to first order. Lemma 4.1. Under Assumption 1, let Z ∈ H. There exist a neighborhood B of Z and Cd > 0 such that ∥d(X) − r(X)∥ ≤ Cd ∥r(X)∥2 for all X ∈ B.
Proof. Let Xp = PH (X) such that r(X) = Xp −X and h(Xp ) = 0. The projection gives X − Xp ∈ NH (Xp ) = im(Dh∗Xp ), and thus −r(X) = Dh∗Xp (v) for some v ∈ Rq . A Taylor expansion yields h(X) = h(Xp ) + DhXp (X − Xp ) + O(∥r(X)∥2 ) = DhXp Dh∗Xp (v) + O(∥r(X)∥2 ). Since DhX = DhXp + O(∥r(X)∥), we obtain d(X) = −Dh∗X (DhX Dh∗X )−1 h(X) = −Dh∗Xp (v) + O(∥r(X)∥2 ) = r(X) + O(∥r(X)∥2 ). Following (3.2), the feasibility direction Gh (X) is given as follows,
(4.2)
Gh (X) = PTM (X) (d(X)),
the projected Gauss–Newton direction. An ensuing question is whether the projected component of d(X) still pulls the iterates towards H. To answer this, we investigate the intersection geometry of H and M, taking into account the intrinsic transversality. Intuitively, to guarantee that the projected Gh (X) inherits the ability of d(X) to decrease dist(·, H), the original d(X) should retain sufficient norm along the tangent space of M, motivating us to propose the following condition.
Definition 4.2 (Projection-based transversality). Two manifolds X , Y ⊆ E satisfy projection-based transversality at Z ∈ X ∩Y if there exist a constant κ′ ∈ (0, 1] and a neighborhood B of Z such that, denoting u = (X − Y )/∥X − Y ∥, ′ (i) for all Y ∈ (Y ∩ B) \ X with X = PX (Y ), one has dist u, NY (Y ) ≥ κ; ′ (ii) for all X ∈ (X ∩ B) \ Y with Y = PY (X), one has dist u, −NX (X) ≥ κ .
Since the above condition constrains only the projection pairs (X, PY (X)) and (PX (Y ), Y ), it is a priori weaker than intrinsic transversality, which constrains all pairs around the common point (see Definition 2.2). Nevertheless, we then show in Theorem 4.4 that, for C 2 submanifolds, the two conditions are equivalent, before which the following lemma is presented as a preliminary. Lemma 4.3. For all nonzero a, b ∈ E, a/∥a∥ − b/∥b∥ ≤ 2∥a − b∥/∥a∥.
Proof. It is seen from the equality a/∥a∥−b/∥b∥ = (a−b)/∥a∥+b(1/∥a∥−1/∥b∥).
Theorem 4.4. Let X and Y be C 2 embedded submanifolds of E, and let Z ∈ X ∩ Y. Then X and Y are intrinsically transversal at Z if and only if they satisfy the projection-based transversality at Z. Proof. For the “only if” part, assume that the intrinsic transversality holds with constant κ. Let Y ∈ (Y ∩ B) \ X and X = PX (Y ), u = (X − Y )/∥X − Y ∥. If X ∈ / Y, the optimality of the projection gives Y − X ∈ NX (X), and thus u ∈ −NX (X), which
10
Y. YANG, B. GAO, AND Y.-X. YUAN
together with (2.5) directly gives dist(u, NY (Y )) ≥ κ. If X ∈ X ∩ Y, then X, Y ∈ Y and Lemma 3.2 applied with base point Y gives dist(X − Y, TY (Y )) ≤ CY ∥X − Y ∥2 for some CY > 0. Hence dist(u, TY (Y )) ≤ CY ∥X − Y ∥ and dist(u, NY (Y )) ≥ 1/2 for ∥X − Y ∥ small enough. Consequently, condition (i) holds with κ′ = min{κ, 1/2}. The symmetric argument yields (ii). For the “if” part, assume that projection-based transversality holds with constant κ′ > 0. By Lemma 3.3, for S ∈ {X , Y}, there exist a neighborhood B of Z, δ > 0, and L > 0 such that PS is single-valued and L-Lipschitz on {X̃ ∈ E | dist(X̃, S ∩ B) < δ}, and PS (X + tv) = X for X ∈ S ∩ B, unit vector v ∈ NS (X), and t ∈ (0, δ). Suppose, for contradiction, that intrinsic transversality fails at Z. Then for every κ ∈ (0, 1), there exist X ∈ (X ∩ B) \ Y and Y ∈ (Y ∩ B) \ X with (4.3)
dist(u, NY (Y )) < κ and
dist(u, −NX (X)) < κ,
where u = (X − Y )/∥X − Y ∥. We show that (4.3) indeed implies a violation of projection-based transversality for the projection pair (X, Y ⋆ ) with Y ⋆ := PY (X). Set ε := dist(u, NY (Y )), τ := ∥X − Y ∥, and w := PNY (Y ) (u). Since ε < 1, the vector w is nonzero, and thus v := w/∥w∥ is a well-defined unit√vector in NY (Y ). 2 The orthogonal √ decomposition u = w + (u − w) yields ∥w∥ = 1 − ε and ∥u − v∥2 = 2 − 2 1 − ε2 ≤ 2ε2 . Shrinking B if necessary such that τ < δ, Lemma 3.3 ensures PY (Y + τ v) = Y . Since X = Y + τ u, the Lipschitz √ continuity of PY yields ∥Y ⋆ − Y ∥ = ∥PY (X) − PY (Y + τ v)∥ √ ≤ Lτ ∥u − v∥ ≤ 2 Lτ ε, and consequently, ∥X − Y ⋆ ∥ ≥ τ − ∥Y − Y ⋆ ∥ ≥ τ (1 − 2Lε) ≥ τ /2 once ε is small enough. Let u⋆ := (X − Y ⋆ )/∥X − Y ⋆ ∥. Lemma 4.3 with a = X − Y ⋆ and b = X − Y gives √ ⋆ ⋆ ∥u − u∥ ≤ 4∥Y − Y ∥/τ ≤ 4 2Lε. The 1-Lipschitz property of the distance √ to a ⋆ ⋆ subspace then gives dist(u , −N (X)) ≤ dist(u, −N (X)) + ∥u − u∥ < κ + 4 2Lε < X X √ (1 + 4√ 2L)κ, where the inequalities follow from (4.3). When κ is small such that (1 + 4 2L)κ < κ′ , it contradicts projection-based transversality for (X, Y ⋆ ). Applying Theorem 4.4 to (X , Y) = (M, H), intrinsic transversality ensures that the projection direction from M onto H is not nearly normal to M, that is, ∥PTM (X) (r(X)/∥r(X)∥) ∥ = dist (r(X)/∥r(X)∥, NM (X)) ≥ κ′ . Recalling that the Gauss–Newton direction d(X) approximates the projection direction r(X), a similar property carries over to Gh (X) = PTM (X) (d(X)).
Proposition 4.5. Assume that H and M are intrinsically transversal at Z ∈ H ∩ M with constant κ. Then there exist a neighborhood B of Z and a constant κ0 > 0 such that for all X ∈ M ∩ B, dist (d(X), NM (X)) ≥ κ0 ∥d(X)∥. Equivalently, (4.4)
∥Gh (X)∥ = ∥PTM (X) (d(X))∥ ≥ κ0 ∥d(X)∥.
Proof. By Lemma 3.3, PH is single-valued and Lipschitz in a neighborhood B of Z. Take X ∈ M ∩ B \ H and set Xp := PH (X), r(X) := Xp − X, r̂ := r(X)/∥r(X)∥. The projection gives X − Xp ∈ NH (Xp ), which implies that −r̂ ∈ NH (Xp ). We first show that r̂ is bounded away from NM (X). If Xp ∈ / M, applying intrinsic transversality (2.5) to the pair (X, Xp ) with u = −r̂ yields dist(−r̂, −NM (X)) ≥ κ, i.e., dist(r̂, NM (X)) ≥ κ. If Xp ∈ M, then both X and Xp lie on M, and Lemma 3.2, after shrinking B if necessary, gives dist(r̂, NM (X)) ≥ 1/2. In either case, we have (4.5)
dist(r̂, NM (X)) ≥ min{κ, 1/2} =: κ1 .
OPTIMIZATION OVER THE INTERSECTION OF MANIFOLDS
11
It remains to transfer this bound from the projection residual r(X) to the Gauss– Newton direction d(X). By Lemma 4.1, ∥d(X) − r(X)∥ ≤ C∥r(X)∥2 for a constant C > 0. Shrinking B such that ∥d(X) − r(X)∥ ≤ 21 ∥r(X)∥ ensures that ∥d(X)∥ ≥ 1 ˆ ˆ 2 ∥r(X)∥ > 0. Writing d := d(X)/∥d(X)∥, Lemma 4.3 gives ∥d − r̂∥ ≤ 2∥d(X) − r(X)∥/∥d(X)∥ ≤ 4C∥r(X)∥. This, together with (4.5) and the 1-Lipschitz property ˆ NM (X)) ≥ κ1 − 4C∥r(X)∥. Shrinking of the distance to a subspace, yields dist(d, ˆ NM (X)) ≥ κ1 /2. The B further such that 4C∥r(X)∥ < κ1 /2, we obtain dist(d, conclusion follows by taking κ0 := κ1 /2 and noting ∥PTM (X) (d)∥ = dist(d, NM (X)).
Proposition 4.5 shows that the projected Gauss–Newton direction preserves a certain component of the original d(X), and thus Gh (X) serves as a descent direction for the feasibility measure dist(·, H); the formal statement is deferred to Proposition 6.5. 5. Clean intersection and optimality direction. We now turn to the optimality direction. Recall from (3.2) that constructing Gf reduces to choosing η(X) ∈ E and projecting it onto S(X) = THX ∩M (X). A natural choice for decreasing the objective is the negative gradient η(X) = −∇f (X), yielding (5.1)
Gf (X) = PS(X) (−∇f (X)).
5.1. Equivalence between intrinsic transversality and clean intersection. The immediate obstacle is that the structure of the tangent cone THX ∩M (X) remains unclear, which impedes the computation of (5.1). To circumvent it, we then investigate the relationship between THX ∩M (X) and the intersection THX (X) ∩ TM (X), enlightened by the general inclusion (2.2). As a result, it is concluded in the following theorem that intrinsic transversality indeed implies clean intersection (see Definition 2.1) for smooth manifolds, answering the question posed in [21, §8]. Theorem 5.1. Let X , Y ⊆ E be C 1 embedded submanifolds and let Z ∈ X ∩ Y. Then X and Y are intrinsically transversal at Z if and only if they intersect cleanly at Z, i.e., there exists a neighborhood B of Z such that X ∩ Y ∩ B is a C 1 embedded submanifold and (5.2)
TX ∩Y (X) = TX (X) ∩ TY (X)
for all X ∈ X ∩ Y ∩ B.
Proof. The “if” part is shown in [21, §3], and thus it suffices to prove the converse. By an appropriate translation, we assume Z = 0 and write L := TX (0). According to Lemma 3.4, there exist a neighborhood UL ⊆ L of 0 and a C 1 map µ : UL → L⊥ with X ∩ B1 = { u + µ(u) | u ∈ UL }, µ(0) = 0, and Dµ(0) = 0 for a neighborhood B1 of Z. Define the C 1 map Ψ : E → E by Ψ(u + n) := u + (n − µ(u)) for u ∈ L, n ∈ L⊥ . One readily verifies that Ψ(0) = 0 and DΨ(0) = I, and thus Ψ is a local C 1 diffeomorphism near 0. Moreover, Ψ(u + µ(u)) = u for all u ∈ UL , and thus Ψ maps X locally onto L. By Lemma 3.6, L and Ψ(Y) are intrinsically transversal at 0. It therefore suffices to prove clean intersection for L ∩ Ψ(Y). For brevity, we slightly abuse the notation by renaming Ψ(Y) as Y in the following analysis. Set Q := TY (0) and decompose E into four orthogonal subspaces by E = K ⊕ U ⊕ V ⊕ W, where K := L ∩ Q, L = K ⊕ U, Q = K ⊕ V, and W := (L + Q)⊥ ; we remark that the orthogonality between U and V can be guaranteed after applying an invertible linear transformation to (L, Y) if needed. Points are written as p = (k, u, v, w) accordingly. By Lemma 3.4 applied to Y with tangent space Q, there exists a C 1 map g : UQ → Q⊥ = U ⊕ W with g(0) = 0 and Dg(0) = 0 such that, writing g(k, v) = (a(k, v), b(k, v)) with a : K ⊕ V → U and b : K ⊕ V → W, (5.3)
Y ∩ B3 = { (k, a(k, v), v, b(k, v)) | (k, v) ∈ UQ }
12
Y. YANG, B. GAO, AND Y.-X. YUAN
for a neighborhood B3 of 0. Throughout the proof, we identify the tuple (k, u, v, w) ∈ K ⊕ U ⊕ V ⊕ W with the element k + u + v + w ∈ E, according to the decomposition. We claim that there exists ε > 0 such that (5.4)
b(k, 0) = 0
for all k ∈ K with ∥k∥ < ε.
Suppose otherwise: there exists a sequence kj → 0 with b(kj , 0) ̸= 0 for all j. Setting v = 0 in (5.3) produces the pairs xj := (kj , a(kj , 0), 0, 0) ∈ L, yj := (kj , a(kj , 0), 0, b(kj , 0)) ∈ Y. Since b(kj , 0) ̸= 0, the point yj has a nonzero Wcomponent, and thus yj ∈ / L. Moreover, PQ (xj ) = (kj , 0) = PQ (yj ), and since PQ |Y is a local diffeomorphism, yj is the unique point of Y near 0 projecting to (kj , 0); as xj ̸= yj , we have xj ∈ / Y. The difference xj − yj = (0, 0, 0, −b(kj , 0)) lies in W ⊆ NL (xj ) = V ⊕ W, and thus the unit vector uj := (xj − yj )/ ∥xj − yj ∥ satisfies dist(uj , −NL (xj )) = 0. On the other hand, since Y is locally the graph t 7→ t + g(t) over Q, the tangent space at yj is TY (yj ) = {ṫ + Dg(kj , 0)[ṫ] | ṫ ∈ Q}. Consider the vector nj := −Dg(kj , 0)∗ uj + uj , and we can verify that ⟨nj , ṫ + Dg(kj , 0)[ṫ]⟩ = 0 for all ṫ ∈ Q, thereby nj ∈ NY (yj ). Noticing that Dg(kj , 0) → 0, we have ∥nj − uj ∥ → 0, and thus limj→∞ max {dist(uj , NY (yj )), dist(uj , −NL (xj ))} = 0, which contradicts the intrinsic transversality of L ∩ Y at 0. According to (5.3) and (5.4), the intersection L ∩ Y is parametrized by F : K → E : ξ 7→ (ξ, a(ξ, 0), 0, 0) in the sense that L ∩ Y ∩ B4 = {F (ξ) | ξ ∈ K, ∥ξ∥ < ε} for a neighborhood B4 of 0. Note that F is C 1 with the injective differential DF (0) : ξ 7→ (ξ, 0, 0, 0) (since Da(0, 0) = 0). By [32, Theorem 4.25], F is a local embedding, and thus L ∩ Y is a C 1 submanifold of dimension dim(K) near the origin. It remains to verify TL∩Y (p) = TL (p) ∩ TY (p) at every nearby intersection point p = F (ξ). The inclusion “⊆” holds in general by (2.2). For “⊇”, take any (k̇, u̇, v̇, ẇ) ∈ TL (p) ∩ TY (p). Membership in TL (p) = K ⊕ U forces v̇ = 0 and ẇ = 0. Differentiating (5.3), (k̇, u̇, v̇, ẇ) ∈ TY (p) requires u̇ = ∂k a(ξ, 0)[k̇]+∂v a(ξ, 0)[v̇] and ẇ = ∂k b(ξ, 0)[k̇] + ∂v b(ξ, 0)[v̇]. Substituting v̇ = 0 gives u̇ = ∂k a(ξ, 0)[k̇] and ẇ = ∂k b(ξ, 0)[k̇]. Since (5.4) implies ∂k b(ξ, 0) = 0, the condition ẇ = 0 is satisfied, and the vector (k̇, ∂k a(ξ, 0)[k̇], 0, 0) = DF (ξ)[k̇] lies in im(DF (ξ)) = TL∩Y (p), establishing the “⊇”. Applying Lemma 3.5 through Ψ−1 completes the proof. 5.2. Projection onto the tangent space. Recalling from Assumption 2 that in the neighborhood K, the manifolds HX and M are intrinsically transversal at the common point X. Therefore, applying Theorem 5.1 to the pair (HX , M) gives the formula of the tangent cone (1.1). This identification reveals that S(X) = THX ∩M (X) is a linear subspace, determined by TM (X) and THX (X). Subsequently, the projection of any ξ ∈ E onto S(X) admits a closed-form expression; see the following proposition. Proposition 5.2. Assume that HX and M are intrinsically transversal at X and let ΦM (X) := PTM (X) ◦ Dh∗X . Then for all ξ ∈ E, it holds that (5.5)
† PS(X) (ξ) = PTM (X) (ξ) − ΦM (X) DhX ◦ ΦM (X) DhX PTM (X) (ξ) .
If additionally, HX and M are transversal at X, then DhX ◦ ΦM (X) is invertible and the Moore–Penrose pseudoinverse reduces to the ordinary inverse. Proof. We note that THX (X) = ker(DhX ) and NHX (X) = im(Dh∗X ). By Theorem 5.1, the clean intersection (1.1) holds. This implies that S(X) is a linear subspace ¯ with ξ¯ := PT (X) (ξ). Therefore, it suffices of TM (X), and thus PS(X) (ξ) = PS(X) (ξ) M ¯ to project ξ ∈ TM (X) onto S(X).
OPTIMIZATION OVER THE INTERSECTION OF MANIFOLDS
13
We claim that TM (X) = S(X) ⊕ im(ΦM (X)). To see this, denote the restricted linear map L := DhX |TM (X) : TM (X) → Rq with S(X) = ker(L). The adjoint L∗ : Rq → TM (X) is characterized by ⟨Lη, v⟩ = ⟨η, L∗ v⟩ for all η ∈ TM (X) and v ∈ Rq . Since ⟨DhX (η), v⟩ = ⟨η, Dh∗X (v)⟩ = ⟨η, PTM (X) (Dh∗X (v))⟩, where the last equality uses η ∈ TM (X), we identify L∗ (v) = ΦM (X)(v), and thus im(ΦM (X)) = im(L∗ ). The theory of linear algebra indicates that ker(L) and im(L∗ ) are orthogonal comple¯ ments within TM (X). Hence, the projection of ξ¯ onto S(X) is ξ¯ − Pim(ΦM (X)) (ξ). q ¯ To compute Pim(ΦM (X)) (ξ), write it as ΦM (X)λ for some λ ∈ R . Since ξ¯ − ¯ The ΦM (X)λ ∈ S(X) = ker(DhX |TM (X) ), we have DhX (ΦM (X)λ) = DhX (ξ). ¯ right-hand side satisfies DhX (ξ) ∈ im(DhX |TM (X) ) = im(DhX ◦ΦM (X)), and thus the ¯ which system is consistent. A representative solution is λ = (DhX ◦ ΦM (X))† DhX (ξ), ¯ = ξ¯−ΦM (X)λ to obtain (5.5). Under is then substituted into the projection PS(X) (ξ) transversality, TM (X)+ker(DhX ) = E, equivalent to NM (X)∩im(Dh∗X ) = {0}. This means ΦM (X) = PTM (X) ◦ Dh∗X is injective, and thus DhX ◦ ΦM (X) is invertible. Proposition 5.2 provides a recipe for computing the projection onto the tangent space of the manifold intersection under intrinsic transversality, and in particular yields the optimality direction Gf (X) (5.1). We then apply it to the case when H and M are the hyperboloid and the low-rank manifold, respectively; see Appendix A. 6. Geometric method and convergence analysis. Combining the feasibility direction (4.1) and the optimality direction (5.1), we now present the geometric method via orthogonal tangent directions in Algorithm 1. Algorithm 1 Geometric method via Orthogonal Tangent Directions (GOTD) Input: Initial point X0 ∈ M, iteration number K, step sizes {αk , βk }K−1 k=0 1: for k = 0, 1, . . . , K − 1 do −1 ∗ 2: Gauss–Newton direction dk = −Dh∗ h(Xk ) Xk DhXk DhXk 3: Feasibility direction Gh (Xk ) = PTM (Xk ) (dk ) 4: Optimality direction Gf (Xk ) = PSk (−∇f (Xk )), Sk := THX (Xk ) ∩ TM (Xk ) k 5: Update Xk+1 = RM Xk αk Gh (Xk ) + βk Gf (Xk ) 6: end for Output: XK Noting that the iterates generated by Algorithm 1 are not necessarily on the set H, we carry out the convergence analysis in a neighborhood of H and show that the iterates approach H asymptotically. In light of the tubular neighborhood theorem [32, Theorem 6.24], we make the following assumption. Assumption 3 (Tubular neighborhood). There exists ∆ > 0 such that PH is single-valued and Lipschitz on the tube K := {X ∈ E | δ(X) < ∆}, where δ(X) := dist(X, H). In addition, finf := inf{f (X) | X ∈ K ∩ M} > −∞.
We also require some regularity conditions on h, commonly adopted in equalityconstrained optimization [44, 54, 2]. Assumption 4. There exist constants σ, µ, Lh > 0 such that σ 2 Iq ⪯ DhX Dh∗X ⪯ µ Iq and ∥DhX − DhY ∥ ≤ Lh ∥X − Y ∥ for all X, Y ∈ K. 2
Assumption 5. The projection PH is C 1 on K with Lipschitz derivative: there exists Lp > 0 such that ∥D(PH )X − D(PH )Y ∥ ≤ Lp ∥X − Y ∥ for all X, Y ∈ K.
Turning to M, the next assumption takes into account the approximation prop-
14
Y. YANG, B. GAO, AND Y.-X. YUAN
erty of the retraction and the Riemannian smoothness of f , both standard in Riemannian optimization [13]. Assumption 6. There exist LR , Lf , Mf , ρM > 0 such that for all X ∈ K ∩ M and 2 η ∈ TM (X) with ∥η∥ ≤ ρM , ∥X + η − RM X (η)∥ ≤ LR ∥η∥ , ∥gradM f (X)∥ ≤ Mf , and (6.1)
Lf ∥η∥2 . f RM X (η) ≤ f (X) + ⟨gradM f (X), η⟩ + 2
Assumptions 4–5 and Assumption 6 consider the regularities of H and M, respectively. We then turn to the intersection geometry. Specifically, recall that Assumption 2 concerns the intrinsic transversality of HX ∩ M in K; this, in views of Proposition 4.5, justifies the following assumption with a uniform constant κ0 > 0. Assumption 7. There exists κ0 > 0 such that for all X ∈ K ∩ M, (6.2)
∥PTM (X) (d(X))∥ ≥ κ0 ∥d(X)∥.
6.1. Auxiliary lemmas. This section collects properties regarding the level sets defined by h in K, and thus all the lemmas below are based on Assumptions 1, 3, 4, and 5. We write r(X) = PH (X) − X as the projection residual. The following lemma presents the error bound condition of h and shows that the Gauss–Newton direction d(X) (4.1) approximates r(X) to second order. Lemma 6.1. There exist constants ch , Ch , Cd > 0 such that ch δ(X) ≤ ∥h(X)∥ ≤ Ch δ(X) and ∥d(X) − r(X)∥ ≤ Cd δ(X)2 for all X ∈ K.
Proof. Let Y := PH (X), such that h(Y ) = 0. Since h is C 2 , we have h(X) = h(Y ) + DhY (X − Y ) + R where ∥R∥ ≤ L2h ∥X − Y ∥2 = L2h δ(X)2 . Noting that X − Y = −r(X), we have h(X) = −DhY [r(X)] + R. Since the projection gives r(X) ∈ NH (Y ) = im(Dh∗Y ), Assumption 4 applied to DhY |im(Dh∗Y ) yields σδ(X) ≤ ∥DhY [r(X)]∥ ≤ µδ(X), and thus ∥h(X)∥ ≤ µδ(X)+ L2h δ(X)2 and ∥h(X)∥ ≥ σδ(X)− Lh Lh 2 2 δ(X) . Shrinking ∆ such that 2 ∆ ≤ σ/2 yields the inequalities with ch = σ/2 Lh and Ch = µ + 2 ∆. The second estimate follows from the analysis of Lemma 4.1.
The next lemma relates the magnitude of the feasibility direction Gh with the distance to the feasible region H. Lemma 6.2. There exist constants cn , Cn > 0 such that for all X ∈ K,
(6.3)
cn δ(X) ≤ ∥Gh (X)∥ ≤ Cn δ(X).
Proof. Let Φ(X) := Dh∗X (DhX Dh∗X )−1 . By Assumption 4, ∥y∥/µ ≤ ∥Φ(X)y∥ ≤ ∥y∥/σ for all y ∈ Rq . For the lower bound, Assumption 7 gives ∥Gh (X)∥ ≥ κ0 ∥d(X)∥. Since ∥d(X)∥ = ∥Φ(X)h(X)∥ ≥ ∥h(X)∥/µ ≥ (ch /µ)δ(X) by Lemma 6.1, we obtain cn = κ0 ch /µ. For the upper bound, ∥Gh (X)∥ ≤ ∥d(X)∥ = ∥Φ(X)∥∥h(X)∥ ≤ (Ch /σ)δ(X), admitting Cn = Ch /σ. We then consider the derivative of the map r; the first observation is an identity that simplifies the inner product involving Dr. Lemma 6.3. For all X ∈ K and η ∈ E, one has ⟨r(X), DrX [η]⟩ = −⟨r(X), η⟩.
Proof. Since DrX = D(PH )X − I and PH maps into H, we have D(PH )X [η] ∈ TH (PH (X)). The identity then follows from r(X) ∈ NH (PH (X)).
This section ends with a bound on the operator norm of Dr, together with a sharper estimate along directions tangent to H.
OPTIMIZATION OVER THE INTERSECTION OF MANIFOLDS
15
Lemma 6.4. For all X ∈ K with Y := PH (X), ∥DrX ∥ ≤ 1+Lp ∆ and ∥DrX [v]∥ ≤ Lp δ(X) ∥v∥ for all v ∈ TH (Y ). If additionally Assumption 4 holds, then there exists a constant Cr > 0 such that ∥DrX [Gf ]∥ ≤ Cr δ(X) ∥Gf ∥ for all Gf ∈ S(X).
Proof. Since PH |H = Id, D(PH )Y acts as the identity on TH (Y ) and vanishes on NH (Y ), giving ∥DrY ∥ = 1. Assumption 5 then yields ∥DrX ∥ ≤ 1 + Lp δ ≤ 1 + Lp ∆. For v ∈ TH (Y ), DrY [v] = 0, and thus ∥DrX [v]∥ = ∥DrX [v] − DrY [v]∥ ≤ Lp δ∥v∥. 0 For the Cr bound, decompose Gf = G0f + G⊥ f with Gf := Pker(DhY ) (Gf ) and 0 ⊥ G⊥ f := Gf − Gf ∈ ker(DhY ) . The observation DhX [Gf ] = 0 and the Lipschitz property of Dh indicate that ∥DhY [Gf ]∥ ≤ Lh δ∥Gf ∥; as DhY [G⊥ f ] = DhY [Gf ], it 0 follows that ∥G⊥ ∥ ≤ (L /σ)δ∥G ∥. Since G ∈ T (Y ), the tangent bound gives h f H f f 0 ⊥ ⊥ ∥DrX [Gf ]∥ ≤ Lp δ∥Gf ∥. Moreover, note that Gf ∈ ker(DhY ) = NH (Y ). Hence, ⊥ ⊥ ⊥ ⊥ D(PH )Y [G⊥ f ] = 0 and DrY [Gf ] = −Gf . By Assumption 5, ∥DrX [Gf ] − DrY [Gf ]∥ ≤ ⊥ ⊥ Lp δ∥G⊥ f ∥, implying that ∥DrX [Gf ]∥ ≤ (1 + Lp ∆)∥Gf ∥ ≤ (1 + Lp ∆)(Lh /σ)δ∥Gf ∥. The triangle inequality concludes that Cr can be chosen as Lp + (1 + Lp ∆)Lh /σ.
6.2. One-step estimates of feasibility and optimality. We then derive the main ingredients of the analysis: one-step estimates for the feasibility measure δ and the objective value f under the update rule (3.1). We begin with the following proposition, which reveals that the feasibility direction Gh drives the decrease in δ, while the optimality direction Gf introduces only higher-order perturbations. Proposition 6.5. Under Assumptions 1–7, there exist ᾱ, β̄, a1 , a2 , a3 , a4 > 0 such that for all X ∈ K ∩ M, 0 < α ≤ ᾱ, and 0 < β ≤ β̄, the iterate X+ = RM X (α Gh (X) + β Gf (X)) satisfies (6.4)
δ(X+ ) ≤ (1 − a1 α + a2 α2 ) δ(X) + a3 β δ(X) ∥Gf (X)∥ + a4 β 2 ∥Gf (X)∥2 .
Proof. Write η := α Gh + β Gf , δ := δ(X), δ+ := δ(X+ ), Y := PHX , and κ1 := (κ0 ch /(2µ))2 . We restrict ∆ ≤ min{κ0 ch /(2µCd ), κ1 /(4Cd )} and set ᾱ = min{1/(2κ1 ), ρ′ /(2Cn ∆)}, β̄ = ρ′ /(2Mf ) with ρ′ := min{ρM , 1/LR }, so that ∥η∥ ≤ αCn ∆ + βMf ≤ ρ′ and Assumption 6 applies: X+ = X + η + eR with ∥eR ∥ ≤ LR ∥η∥2 . A Taylor expansion of r at X applied to the displacement η + eR gives L δ+ ≤ ∥r(X) + DrX [η + eR ]∥ + 2p ∥η + eR ∥2 . By Lemma 6.4, ∥DrX ∥ ≤ 1 + Lp ∆. In 2 addition, ∥eR ∥ ≤ LR ∥η∥ ≤ ∥η∥ (by ∥η∥ ≤ 1/LR ) implies ∥η + eR ∥ ≤ 2∥η∥, and thus (6.5)
δ+ ≤ ∥r(X) + DrX [η]∥ + cR ∥η∥2 ,
with cR := (1 + Lp ∆)LR + 2Lp .
Expanding ∥r + Dr[η]∥2 according to Lemma 6.3, (6.6)
∥r + Dr[η]∥2 = δ 2 − 2α⟨r, Gh ⟩ − 2β⟨r, Gf ⟩ + ∥Dr[η]∥2 .
We estimate the terms on the right. For the term involving Gh , we note that Gh = PTM (X) (d) reveals that ⟨r, Gh ⟩ = ⟨PTM (X) (r), d⟩, and writing d = r + (d − r) with ∥d − r∥ ≤ Cd δ 2 (Lemma 6.1) yields ⟨r, Gh ⟩ = ∥PTM (X) (r)∥2 + ⟨PTM (X) (r), d − r⟩ ≥ ∥PTM (X) (r)∥2 −Cd δ 3 . The triangle inequality and Assumption 7 give ∥PTM (X) (r)∥ ≥ κ0 ∥d∥−Cd δ 2 ≥ (κ0 ch /µ−Cd δ)δ; restricting ∆ ≤ κ0 ch /(2µCd ) ensures ∥PTM (X) (r)∥ ≥ κ0 ch δ/(2µ), and thus ⟨r, Gh ⟩ ≥ κ1 δ 2 − Cd δ 3 . For the term involving Gf , the orthogonality ⟨d, Gf ⟩ = 0 (Lemma 3.1) gives |⟨r, Gf ⟩| = |⟨r − d, Gf ⟩| ≤ Cd δ 2 ∥Gf ∥. For the last term, Lemma 6.2 and Lemma 6.4 yield ∥Dr[η]∥2 ≤ 2(1 + Lp ∆)2 Cn2 α2 δ 2 + 2Cr2 β 2 δ 2 ∥Gf ∥2 . Substituting the estimates into (6.6) shows ∥r + Dr[η]∥2 ≤ δ 2 (1 − t) with t := 2α(κ1 − Cd δ) − 2βCd ∥Gf ∥ − 2(1 + Lp ∆)2 Cn2 α2 − 2Cr2 β 2 ∥Gf ∥2 . The restriction on ∆ ensures κ1 − Cd δ ≥ κ1 /2, and ᾱ ≤ 1/(2κ1 ) gives 2ακ1 ≤ 1, hence t ≤ 1.
16
Y. YANG, B. GAO, AND Y.-X. YUAN
√ Applying 1 − t ≤ 1 − t/2 leads to ∥r + Dr[η]∥ ≤ δ − κ1 αδ + Cd αδ 2 + Cd βδ∥Gf ∥ + (1 + Lp ∆)2 Cn2 α2 δ + Cr2 β 2 δ∥Gf ∥2 . Adding cR ∥η∥2 ≤ 2cR α2 Cn2 δ 2 + 2cR β 2 ∥Gf ∥2 from (6.5) and using δ ≤ ∆ to replace each δ 2 by ∆δ, we conclude with a1 = 3κ4 1 , a2 = (1 + Lp ∆)2 Cn2 + 2cR Cn2 ∆, a3 = Cd , and a4 = Cr2 ∆ + 2cR . We then resort to the Riemannian smoothness assumption over M, proving that the optimality direction Gf delivers the descent property of the objective f . Proposition 6.6. Under Assumptions 1–7, let X ∈ K ∩ M, 0 < α ≤ ᾱ, 0 < β ≤ β̄ (with ᾱ, β̄ from Proposition 6.5), and X+ = RM X (α Gh (X) + β Gf (X)). Then Lf Lf 2 (6.7) f (X+ ) − f (X) ≤ −β 1 − β ∥Gf (X)∥2 + α Mf ∥Gh (X)∥ + α ∥Gh (X)∥2 . 2 2 Proof. Let η := α Gh +β Gf . The constructed ᾱ and β̄ ensure that ∥η∥ ≤ ρM , and L Assumption 6 produces f (X+ ) ≤ f (X) + ⟨gradM f (X), η⟩ + 2f ∥η∥2 . For the inner product, Gf = PS(X) (−gradM f (X)) ∈ TM (X) yields ⟨gradM f, Gf ⟩ = −∥Gf ∥2 , while |⟨gradM f, Gh ⟩| ≤ Mf ∥Gh ∥. The orthogonality ⟨Gh , Gf ⟩ = 0 (Lemma 3.1) gives ∥η∥2 = α2 ∥Gh ∥2 + β 2 ∥Gf ∥2 . Substitute the terms and rearrange them to yield (6.7).
Drawing on Proposition 6.5, it is shown that the iterates are confined to the tube K under appropriate step sizes.
Lemma 6.7. Under Assumptions 1–7, let 0 < ∆0 < ∆. There exist constants ᾱ′ , β̄ , τ > 0 such that for all X ∈ K ∩ M with δ(X) ≤ ∆0 , 0 < α ≤ ᾱ′ , 0 < β ≤ β̄ ′ , and β ≤ τ α, the iterate X+ = RM X (αGh + βGf ) satisfies δ(X+ ) ≤ ∆0 . ′
Proof. Let ᾱ, β̄, a1 , . . . , a4 be as in Proposition 6.5. Using ∥Gf ∥ ≤ Mf and β ≤ τ α in (6.4) gives a3 βδ∥Gf ∥ ≤ a3 τ αδMf and a4 β 2 ∥Gf ∥2 ≤ a4 τ 2 α2 Mf2 . Set τ = a1 /(4a3 Mf ) so that a3 τ Mf = a1 /4, and define ᾱ′ = min ᾱ, a1 /(4a2 ), a1 ∆0 /(2a4 τ 2 Mf2 ) , β̄ ′ = min{β̄, τ ᾱ′ }. Then δ(X+ ) ≤ (1 − a1 α/2)δ + a4 τ 2 α2 Mf2 . For δ ≤ ∆0 , the right-hand side is at most (1 − a1 α/2)∆0 + a4 τ 2 α2 Mf2 ≤ ∆0 , where the last step uses a4 τ 2 αMf2 ≤ a1 ∆0 /2.
By induction, if X0 ∈ M with δ(X0 ) ≤ ∆0 and the step sizes satisfy 0 < αk ≤ ᾱ′ , 0 < βk ≤ β̄ ′ , βk ≤ τ αk for all k, then the iterates generated by Algorithm 1 remain in the tube K, i.e., δ(Xk ) ≤ ∆0 < ∆ for all k ≥ 0. 6.3. Lyapunov function and iteration complexity. We consider the following Lyapunov function to treat the feasibility and the optimality measures, Lλ (X) := f (X) + λ δ(X), with λ > 0 as a balance factor.
(6.8)
Combining the estimates in Propositions 6.5 and 6.6 points to the complexity analysis. Theorem 6.8. Under Assumptions 1–7, let 0 < ∆0 < ∆ and X0 ∈ M with δ(X0 ) ≤ ∆0 . Let λ = 2Mf Cn /a1 . There exist ᾱ′′ , β̄ ′′ , τ ′′ , cδ , ct > 0 such that for constant step sizes α ≤ ᾱ′′ , β ≤ β̄ ′′ , β ≤ τ ′′ α, the iterates of Algorithm 1 satisfy Xk ∈ K ∩ M and Lλ (Xk+1 ) ≤ Lλ (Xk ) − cδ α δ(Xk ) − ct β ∥Gf (Xk )∥2 , for all k ≥ 0. Moreover, for all K ≥ 1, we have (6.9)
min
0≤k≤K−1
δ(Xk ) ≤
Lλ (X0 ) − finf , cδ αK
and δ(Xk ) → 0, ∥Gf (Xk )∥ → 0 as k → ∞.
min
0≤k≤K−1
∥Gf (Xk )∥2 ≤
Lλ (X0 ) − finf , ct βK
OPTIMIZATION OVER THE INTERSECTION OF MANIFOLDS
17
L
Proof. Set ᾱ′′ = min{ᾱ′ , Mf Cn /(4(λa2 + 2f Cn2 ∆))}, β̄ ′′ = min{β̄ ′ , 1/(2Lf + 4λa4 )}, and τ ′′ = min{τ, Mf Cn /(4λ2 a23 ∆)}. Lemma 6.7 reveals that δ(Xk ) ≤ ∆0 for all k. Fix k and write X := Xk , X+ := Xk+1 , δ := δ(X). Substituting ∥Gh ∥ ≤ L Cn δ (Lemma 6.2) into Proposition 6.6 gives f (X+ ) − f (X) ≤ −β 1 − 2f β ∥Gf ∥2 + L α Mf Cn δ + 2f α2 Cn2 δ 2 . Adding λ times the feasibility estimate δ(X+ ) − δ ≤ −a1 αδ + 2 a2 α δ + a3 βδ∥Gf ∥ + a4 β 2 ∥Gf ∥2 from (6.4) and relaxing δ 2 by ∆δ, we have L Lλ (X+ ) − Lλ (X) ≤ − β 1 − 2f β − λa4 β ∥Gf ∥2 + λ a3 β δ ∥Gf ∥ (6.10) L − λ a1 − Mf Cn α δ + λa2 + 2f Cn2 ∆ α2 δ. Applying Young’s inequality λa3 βδ∥Gf ∥ ≤ β4 ∥Gf ∥2 + λ2 a23 ∆βδ and β ≤ τ ′′ α lead to L Lλ (X+ ) − Lλ (X) ≤ − β 43 − 2f β − λa4 β ∥Gf ∥2 (6.11) L − λa1 − Mf Cn − λ2 a23 ∆τ ′′ αδ + λa2 + 2f Cn2 ∆ α2 δ. For the ∥Gf ∥2 coefficient, β̄ ′′ ≤ 1/(2Lf + 4λa4 ) ensures Lf β/2 + λa4 β ≤ 1/4, and L thus 34 − 2f β − λa4 β ≥ 12 . For the αδ coefficient, λ = 2Mf Cn /a1 gives λa1 − Mf Cn = Mf Cn , and τ ′′ ≤ Mf Cn /(4λ2 a23 ∆) gives λ2 a23 ∆τ ′′ ≤ Mf Cn /4, hence λa1 − Mf Cn − λ2 a23 ∆τ ′′ ≥ 43 Mf Cn . For the α2 δ term, the choice of ᾱ′′ ensures (λa2 + Lf 2 2 Cn ∆)α ≤ Mf Cn /4. Assembling the estimates, we obtain Lλ (X+ ) − Lλ (X) ≤ −cδ αδ − ct β∥Gf ∥2 with cδ = Mf Cn /2 and ct = 1/2. Telescoping from k = 0 to K − 1 PK−1 leads to k=0 cδ α δ(Xk ) + ct β∥Gf (Xk )∥2 ≤ Lλ (X0 ) − finf < ∞, from which (6.9) and δ(Xk ) → 0, ∥Gf (Xk )∥ → 0 follow. We present the first-order stationarity condition of problem (P) as follows, which is necessary for the local optimality [51, Theorem 6.12].
Definition 6.9. A point X̄ ∈ H ∩ M is called stationary for problem (P) if ⟨∇f (X̄), η⟩ ≥ 0 for all η ∈ TH∩M (X̄), i.e., −∇f (X̄) ∈ NH∩M (X̄), or equivalently, the projected negative gradient vanishes, i.e., PTH∩M (X̄) (−∇f (X̄)) = 0. Therefore, the quantity ∥Gf (X)∥ = ∥PTHX ∩M (X) (−∇f (X))∥ serves as a suitable evaluation of the stationarity. In this view, (6.9) in Theorem 6.8 indeed delivers an √ O(1/ K) convergence rate for both the feasibility and the optimality measures. We conclude by examining the first-order stationarity of the accumulation points.
Corollary 6.10. Under conditions of Theorem 6.8, every accumulation point X̄ ∈ M of {Xk } satisfies h(X̄) = 0. If, additionally, DhX |TM (X) has constant rank near X̄ on M, then X̄ is a first-order stationary point of (P). Proof. By Theorem 6.8, δ(Xk ) → 0 and ∥Gf (Xk )∥ → 0. The continuity of δ and Xkj → X̄ give δ(X̄) = limj δ(Xkj ) = 0, indicating that X̄ ∈ H. For stationarity, denote LX := DhX |TM (X) : TM (X) → Rq . As shown in the proof of Proposition 5.2, the adjoint L∗X : Rq → TM (X) coincides with ΦM (X) = PTM (X) ◦ Dh∗X , thereby DhX ◦ ΦM (X) = LX L∗X ∈ Rq×q . Since rank(LX L∗X ) = rank(LX ), the constant rank condition on LX is equivalent to that on LX L∗X , and thus the pseudoinverse of LX L∗X is continuous. By the formula (5.5) in Proposition 5.2, the map X 7→ PS(X) is then continuous near X̄. Passing to the limit along Xkj → X̄ in Gf (Xkj ) = PS(Xkj ) (−∇f (Xkj )) → 0 gives PS(X̄) (−∇f (X̄)) = 0. Remark 6.11. The constant rank condition on DhX |TM (X) coincides with the relaxed constant rank constraint qualification (RCRCQ) introduced in [6] for Riemannian optimization problems with equality constraints. This, in views of Theorem 5.1,
18
Y. YANG, B. GAO, AND Y.-X. YUAN
amounts to requiring that the slices HX ∩ M have constant dimension as X varies near X̄ on M. Such a property is satisfied by all the instances outlined in Table 1. 7. Numerical experiments. In this section, we evaluate the proposed GOTD (Algorithm 1) on three applications within the scope of formulation (P), each constrained to an intersection of two manifolds. We adopt the Riemannian augmented Lagrangian method (RALM) [38, Algorithm 1] as a baseline, which can handle Riemannian problems with equality and inequality constraints; we run RALM with the publicly available implementation.1 In addition, each experiment also incorporates some other task-specific methods for comparison, where we note that the Riemannian trust-region method on manifolds is run with Manopt’s default settings. The experiments are produced on a workstation that consists of two Intel(R) Xeon(R) Gold 6330 CPUs (at 2.00GHz×28, 42M Cache), 512GB RAM. All the experiments are carried out in MATLAB (Release 9.7.0) on the CPUs, drawing on the Manopt toolbox [14]. The codes of the proposed method are available at https://github.com/UCAS-YanYang. 7.1. Low-rank approximation of spherical data. Finding a low-rank approximation of normalized data plays a crucial role in various applications. Given A ∈ Ob(m, n) := {X ∈ Rm×n | diag(XX ⊤ ) − 1 = 0}, where rows encode data points with unit length, Chu et al. [20] formulated the approximation task as follows, (7.1)
min
X∈Rm×n
s. t.
1 2 ∥PΩ (X − A)∥ 2 X ∈ Ob(m, n) ∩ Mr ,
where Ω ⊆ {1, 2, . . . , m} × {1, 2, . . . , n} represents observed entries and PΩ defines the sampling operator: PΩ (X)(i, j) = X(i, j) if (i, j) ∈ Ω, otherwise PΩ (X)(i, j) = 0. Apart from RALM, we also take into account the method proposed in [65], which parameterizes Ob(m, n) ∩ Mr by a smooth manifold denoted by Mh . Subsequently, Manopt’s Riemannian gradient descent and Riemannian trust-region methods are invoked on Mh for comparison, denoted by Mh -RGD and Mh -RTR, respectively. Following the test of [65], we generate a ground truth A = POb(m,n) (U ∗ Σ ∗ (V ∗ )⊤ ), where U ∗ ∈ St(m, r) and V ∗ ∈ St(n, r) are obtained by sampling entries from the standard normal distribution N (0, 1) and taking the Q-factors of the QR factorizations, and Σ ∗ ∈ Rr×r is a diagonal matrix with entries sampled from the uniform distribution on (0, 1). The oversampling factor is defined by OS := |Ω|/(r(m + n − r)). All the methods share a common initial point X0 = H0 V0⊤ , where V0 is generated in the same way as V ∗ and H0 = POb(m,r∗ ) (H̃0 ), where entries of H̃0 are sampled from N (0, 1). The termination rules are set for respective methods: max{∥Gh (Xk )∥, ∥Gf (Xk )∥} ≤ 10−10 for GOTD, subproblem accuracy of RALM achieves the tolerance 10−8 (see [38, Algorithm 1]), the norm of Riemannian gradient on Mh achieves 10−10 for Mh -RGD and 10−13 for Mh -RTR. Reconstruction quality is measured by the relative test error ∥PΓ (X − A)∥F /∥PΓ (A)∥F on an independent test set Γ with |Γ | = |Ω|. All reported numbers are averaged over five random seeds. We test with the dimension (m, n) = (5000, 6000), the oversampling factors OS ∈ {6, 7, 8, 9, 10}, and the rank parameters r ∈ {8, 9, 10}. For GOTD, we adopt constant step sizes throughout: α = 1 in all experiments and β tuned over the grid {1, 5, 10, 20, 30, 40, 50}. The same tuning strategy applies to the constant step size in Mh -RGD. Table 2 reports the performance of GOTD across different (OS, r), with the true data matrix successfully recovered in every configuration. 1 https://github.com/losangle/Optimization-on-manifolds-with-extra-constraints
19
OPTIMIZATION OVER THE INTERSECTION OF MANIFOLDS
Table 2 Performance of GOTD on spherical data fitting problem with different oversampling factors and rank parameters, each entry averaged over five seeds. r=8
OS
6 7 8 9 10
r=9
r = 10
Test err.
Time
Test err.
Time
Test err.
Time
1.51e−11 9.21e−12 8.80e−12 1.69e−11 2.83e−12
5.41 2.71 4.45 5.69 2.32
9.38e−12 8.05e−12 1.53e−12 6.63e−13 6.84e−13
3.17 2.72 2.83 13.58 13.89
8.40e−12 8.79e−12 6.76e−13 6.54e−12 6.97e−13
5.50 7.06 13.56 12.15 9.96
Moreover, Figure 4 compares the four methods on the scenario (OS, r) = (6, 10), where the feasibility is measured by ∥h(X)∥ = ∥ diag(XX ⊤ ) − 1∥. RALM and Mh RTR exhibit better iteration complexity, since the former solves a Lagrangian subproblem at each update and the latter exploits Hessian information. In terms of running time, GOTD is more efficient, as each of its iterations invokes only first-order oracles and one retraction on Mr , without additional inner subproblems to solve. GOTD
10-10
0
100
200
300
Iteration
400
500
10-5
10
kh(X)k
Test error
Test error
-5
Mh -RTR 100
100
100
10
Mh -RGD
RALM
-5
10-10
10-10
0
5
10
Time (s)
15
20
10-15
0
5
10
15
20
Time (s)
Fig. 4. Low-rank approximation of spherical data with (m, n) = (5000, 6000), r = r∗ = 10, and OS = 6. Left: test error versus iteration. Middle: test error versus wall-clock time. Right: feasibility violation ∥ diag(Xk Xk⊤ ) − 1∥ versus wall-clock time.
7.2. Low-rank approximation of hyperbolic embeddings. Hyperbolic embeddings have shown effectiveness in representing hierarchical data, with applications ranging from graph embeddings to natural language processing [42, 29]. To capture the hyperbolic geometry, the hyperboloid model is introduced as Hn := {y ∈ Rn+1 | ⟨y, y⟩J + 1 = 0, y1 > 0}, where y1 is the first element of the vector, J := Diag(−1, 1, . . . , 1) is the Lorentz signature, and ⟨x, y⟩J := x⊤ Jy denotes the Lorentzian inner product. As the upper sheet of an n-dimensional hyperboloid in Rn+1 , Hn is a smooth manifold with the distance dH (x, y) := arccosh(−⟨x, y⟩J ). Given m data points which have been embedded in Hn as x̄1 , . . . , x̄m , Jawanpuria et al. [29] proposed to seek x1 , . . . , xm ∈ Hn as the approximations, which share a latent low-dimensional subspace—that is, the matrix X = [x1 , . . . , xm ] is low-rank. Stacking the constraints xi ∈ Hn translates to X lying on the matrix hyperboloid, (7.2)
(n+1)×m Hm | diag(Y ⊤ JY ) + 1 = 0, Y1: > 0}, n := {Y ∈ R
where Y1: extracts the first row of Y and the “>” is understood component-wise. Letting r ≪ min{m, n} denote the rank parameter, finding a low-rank approximation
20
Y. YANG, B. GAO, AND Y.-X. YUAN
Table 3 Results of the four methods on low-rank hyperbolic embeddings on the WordNet mammals subtree. The f0 and f are the initial and returned objective value, respectively. For reference, the mean average precision computed through the original embedding x̄1 , . . . , x̄m is 0.9385. Rank
r=5
r = 10
r = 20
Time (s)
∥h(X)∥
Precision
0.434
80.8
6.81e−8
0.8869
0.994
699.2
1.82e−6
0.6637
Mp -RGD
0.439
92.9
1.89e−12
0.8881
Mp -RTR
0.435
78.3
4.14e−12
0.8883
GOTD
0.394
65.9
5.01e−8
0.9039
RALM
0.949
712.9
3.19e−6
0.8117
Mp -RGD
0.402
99.7
3.38e−12
0.8961
Mp -RTR
0.396
82.5
4.14e−12
0.9024
GOTD
0.369
98.8
1.64e−7
0.9003
RALM
0.877
653.9
3.55e−6
0.8463
Mp -RGD
0.380
100.3
3.94e−12
0.8989
Mp -RTR
0.376
73.0
4.89e−12
0.8981
Algorithm
f /f0
GOTD RALM
of hyperbolic embeddings can be formulated as follows,
(7.3)
m X
arccosh − ⟨xi , x̄i ⟩J
min
f (X) =
s. t.
X ∈ Hm n ∩ Mr+1
X∈R(n+1)×m
i=1
2
where xi denotes the i-th column of X, and the objective is the sum of squared distances between xi and x̄i for i = 1, . . . , m. The existing approach for (7.3), developed by Jawanpuria et al. [29], considers the product manifold Mp := St(n, r) × Hm r , together with a smooth map defined by 1 0⊤ (n+1)×m r ϕ : Mp → R : ϕ(U, Z) = Z. 0n U Since ϕ(Mp ) ⊊ Hm n ∩ Mr+1 , it produces a surrogate problem minZ∈Mp f (ϕ(Z)), thereby allowing the employment of Riemannian gradient descent (RGD) and Riemannian trust-region (RTR) methods on Mp . Our method GOTD, in turn, tackles (n+1)×m the coupled constraint X ∈ Hm ; the compun ∩ Mr+1 in the original space R tation of Gf (X) is referred to Appendix A. We also include RALM for comparison, ⊤ interpreting X ∈ Hm n as the equality constraint h(X) = diag(X JX) + 1 = 0. We evaluate the methods on the mammal subtree of the WordNet lexical database [40], with (n, m) = (300, 1170) and rank parameters r ∈ {5, 10, 20}. The nodes are nouns and the edges encode the “is-a” relationship; for instance, the edge between “squirrel” and “rodent” indicates that a squirrel is a rodent. The original embeddings x̄1 , . . . , x̄m are obtained from the implementation [42]. The quality of the approximations is measured by the mean average precision: given the true edge set E, for each pair (u, v) ∈ E, distance dH (xu , xv ) is ranked among the distances {dH (xu , xw ) | (u, w) ∈ / E}, and the resulting precision over all true edges is averaged. Let X̄ = [x̄1 , . . . , x̄m ], and collect the top-r left singular vectors of the submatrix X̄ ′ := X̄(2:n+1, :) to obtain Ur ∈ St(n, r). We denote Z̄ ′ := Ur⊤ X̄ ′ and write zi′
OPTIMIZATION OVER THE INTERSECTION OF MANIFOLDS
21
p for its i-th column. Augmenting each column as z̄i := [ 1 + ∥zi′ ∥2 , zi′⊤ ]⊤ ∈ Hr for i = 1, . . . , m and assembling them into Z̄ = [z̄1 , . . . , z̄m ], we initialize Mp -RGD and Mp -RTR at (Ur , Z̄), and initialize GOTD and RALM at X0 = ϕ(Ur , Z̄). After conducting initial tests, we use α = 1 and β tuned over {0.1, 0.15, 0.2, 0.5, 0.75, 1} for GOTD; the same tuning strategy applies to the constant step size in Mp -RGD. Table 3 shows that the feasibility measure ∥h∥ reaches the order of 10−7 for GOTD, which is acceptable, and stays at machine precision for Mp -RGD and Mp RTR, since the methods are feasible in the sense that ϕ(Mp ) ⊂ Hm n . For a fair comparison, we post-process the output of GOTD by projecting it onto Hm n ∩ Mr+1 via alternating projections, yielding Xgotd , on which we report the “f /f0 ” and the precision. Across the three ranks, GOTD attains the lowest cost f , together with the highest mean average precision at r = 10 and r = 20. A plausible reason for the enhanced performance is the observation Xgotd ∈ / ϕ(Mp ), which indicates that GOTD searches over a larger set. Moreover, the running time of GOTD is competitive with the Riemannian algorithms conducted on Mp . 7.3. Compressed modes in physics. In quantum mechanics, spatially localized solutions to Schrödinger’s equation have attracted interest recently [48]. Such solutions, known as compressed modes, translate into the sparsity of the discretization matrix. Specifically, with X ∈ St(n, p) collecting p orthonormal wave functions, the task of finding compressed modes in physics can be formulated as follows, (7.4)
min
tr(X ⊤ AX)
s. t.
X ∈ St(n, p) ∩ Cs
X∈Rn×p
where Cs = {Y ∈ Rn×p | ∥Y ∥0 = s}, A ∈ Rn×n is a symmetric matrix discretizing the Hamiltonian, and s > 0 is a parameter depicting the sparsity. Hence (7.4) is an instance of (P) with (H, M) = (St(n, p), Cs ), to which GOTD can be applied. The existing approaches turn to the relaxation by adding an ℓ1 regularization of the wave functions [48, 19], yielding (7.5)
min
tr(X ⊤ AX) + µ∥X∥1
s. t.
X ∈ St(n, p).
X∈Rn×p
In the comparison, we implement GOTD and RALM for (7.4), where RALM treats the constraint X ∈ St(n, p) as the equality constraint h(X) = X ⊤ X − Ip = 0. Moreover, the proximal gradient methods ManPG [19, Algorithm 2] and SLPG [39, Algorithm 3] are applied to (7.5), which are able to tackle nonsmooth objectives. Following the setting of [19], we take the 1D free-electron Hamiltonian − 21 ∂x2 on the interval [0, L] with L = 50, discretized on a uniform grid of n = 256 interior points, and seek p = 15 orthonormal wave functions. The relaxation parameter for the two ℓ1 penalty methods is set to µ = 1/30, taken from the publicly available implementation of ManPG.2 Sparsity is measured by the ratio of the non-zero entry, ρ(X) := (np − ∥X∥0 )/(np). To report the sparsity of X, we treat an entry as zero whenever |Xij | < 10−8 for ManPG and SLPG; the iterates of GOTD and RALM are exactly on Cs at every step, and thus the threshold has no bearing on either method. The initial points for ManPG and SLPG are obtained by running a Riemannian sub-gradient descent on problem (7.5), while GOTD and RALM additionally project the points onto Cs 2 https://github.com/chenshixiang/ManPG/tree/master
22
Y. YANG, B. GAO, AND Y.-X. YUAN
as the warm start. The step size for ManPG and SLPG is 1/(2λmax (A)) = L2 /4n2 with λmax (A) denoting the largest eigenvalue of A. Similarly, we take α = 1 and β = L2 /4n2 for GOTD. The methods ManPG and SLPG are terminated when the iterate satisfies |F (Xk ) − F (Xk+1 )| ≤ 10−7 with F (X) := tr(X ⊤ AX) + µ∥X∥1 , and GOTD is terminated when max{∥Gh (Xk )∥, ∥Gf (Xk )∥} ≤ 10−10 . Figure 5 reports the cost tr(X ⊤ AX), the sparsity ρ, and the feasibility measure ∥h(X)∥ = ∥X ⊤ X −Ip ∥ of the iterates. Specifically, ManPG and SLPG rapidly sparsify the iterates while maintaining the cost value, with the sparsity steady at ρ ≈ 0.6; this motivates us to configure the manifold Cs with s = ρ × np and ρ ∈ {0.6, 0.7} in the formulation (7.4). Consequently, the curves reveal that GOTD attains lower cost value under both levels of sparsity, exhibiting the effectiveness and efficiency of the proposed method. In addition, RALM, also targeting the formulation (7.4), approximately matches the cost value returned by GOTD at the same sparsity, but accrues a slightly higher feasibility violation. GOTD (; = 0:6) 0.8
2.55
0.6
2.5 2.45
RALM (; = 0:6)
0
1
2
Time (s)
3
SLPG
10-5
0.4
10-10
0.2
2.4
ManPG 100
kh(X)k
Sparsity
Cost
GOTD (; = 0:7) 2.6
0
0
1
2
Time (s)
3
10-15
0
1
2
3
Time (s)
Fig. 5. Comparison on finding compressed modes with (n, p) = (256, 15).
8. Conclusions and perspectives. In this work, we propose GOTD, a geometric method for optimization over the intersection of two manifolds H ∩ M, under intrinsic transversality. The main principle is to decompose the update into two orthogonal directions tangent to M: the projected Gauss–Newton direction improving the feasibility and the projected anti-gradient direction decreasing the objective. Importantly, constructions of the two components are promoted by two equivalent characterizations of intrinsic transversality in Theorem 4.4 and Theorem 5.1, respectively. Convergence analysis and numerical experiments validate the effectiveness and efficiency of GOTD. The developed theory suggests several avenues worth exploration. First, given the first-order intersection rule (1.1) enlightening our work, it is reasonable to consider the second-order counterpart T2H∩M (X; η) = T2H (X; η) ∩ T2M (X; η) identified in [64] to devise a second-order extension of the proposed GOTD, where T2 denotes the secondorder tangent set. In addition, adapting Algorithm 1 to accommodate stochastic oracles is of interest for large-scale scientific computation where only noisy estimates of ∇f and h are accessible; some relevant techniques can be found in [54, 2]. Moreover, extending GOTD to handle multiple intersecting constraints X ∈ M ∩ H1 ∩ · · · ∩ HN also appears as a potential direction, which can borrow some ideas from [34]. Acknowledgments. The authors would like to thank P.-A. Absil for his comments and insightful suggestions on this work.
OPTIMIZATION OVER THE INTERSECTION OF MANIFOLDS
23
Appendix A. Projection onto THm (X). Given a rank parameter s with n ∩Ms 1 ≤ s ≤ min{(n + 1), m}, we then derive an efficient procedure for computing the projection onto the space THm (X), which is invoked by the experiment in section 7.2. n ∩Ms Let J = Diag(−1, 1, . . . , 1) ∈ R(n+1)×(n+1) , and write Xi for the i-th column of a ma(n+1)×m trix X. The matrix hyperboloid Hm → n defined in (7.2) is the level set of h : R m ⊤ R with hj (X) = Xj JXj + 1 and the differentials DhX (Z) = 2 diag(X ⊤ JZ) and Dh∗X (λ) = 2JX Diag(λ). Moreover, the normal and tangent spaces are NHm (X) = n m (n+1)×m ⊤ {JX Diag(λ) | λ ∈ R } and THm (X) = {Z ∈ R | Xj JZj = 0, j = 1, . . . , m}, n respectively. Regarding the geometry of Ms , we consider the singular value decomposition X = U ΣV ⊤ with U ∈ St(n + 1, s) and V ∈ St(m, s) to give the characterization TMs (X) = {U U ⊤ Z + ZV V ⊤ − U U ⊤ ZV V ⊤ | Z ∈ R(n+1)×m }, and the orthogonal projection PTMs (X) (Z) = U U ⊤ Z + ZV V ⊤ − U U ⊤ ZV V ⊤ ; see the developments in [58, 47]. Building on the geometry of each manifold, we investigate the intersection. m Proposition A.1. At every point X ∈ Hm n ∩ Ms , the manifolds Hn and Ms (n+1)×m intersect transversally, i.e., THm (X)+TMs (X) = R . Consequently, Hm n ∩Ms n is a smooth submanifold of dimension s(n + m − s + 1) − m.
Proof. It suffices to prove that NHm (X)∩NMs (X) = {0} by definition of transvern sality. To see this, take any JX Diag(λ) ∈ NHm (X), and suppose additionally that it n belongs to NMs (X). Note that for any C ∈ Rm×m , we have XC ∈ TMs (X), and thus JX Diag(λ) ∈ NMs (X) gives ⟨JX Diag(λ), XC⟩ = 0, which forces X ⊤ JX Diag(λ) = 0. Taking the diagonal and recalling diag(X ⊤ JX) = −1 yields λ = 0, thereby NHm (X)∩NMs (X) = {0}. The dimension of the intersection follows from dim(Ms )+ n (n+1)×m dim(Hm ) = s(n + m − s + 1) − m. n ) − dim(R
With Proposition A.1 depicting the transversality, we can deliver the computation of the projection onto THm (X) via Proposition 5.2. To exploit the low-rank n ∩Ms structure X = U ΣV ⊤ , we decompose JX = U P + Q with P := U ⊤ JX ∈ Rs×m and Q := (I − U U ⊤ )JX ∈ R(n+1)×m . Following the notation in Proposition 5.2 that ΦMs (X) = PTMs (X) ◦ Dh∗X , we obtain the following identity, (A.1) ΦMs (X)(λ) = 2PTMs (X) JX Diag(λ) = 2 U P Diag(λ) + Q Diag(λ) V V ⊤ , Letting DP := Diag ∥P1 ∥22 , . . . , ∥Pm ∥22 ∈ Rm×m collect the squared column norms of P , a direct computation then gives, for every λ ∈ Rm , DhX ◦ ΦMs (X)[λ] = 4Aλ with A := DP + (Q⊤ Q) ⊙ (V V ⊤ ). Propositions 5.2 and A.1 guarantee that A is invertible. Substituting (A.1) into (5.5) with ξ ∈ R(n+1)×m and η := PTMs (X) (ξ), we obtain (A.2)
PTHm ∩Ms (X) (ξ) = η − U (P Diag(λ)) − (Q Diag(λ)V )V ⊤ , n
m
where λ ∈ R solves the linear system Aλ = b, with the i-th entry of b given by bi := (JXi )⊤ ηi . Although A is of size m × m, it need not be assembled inPadvance. s m In fact, the action on any w ∈ R admits the factored form Aw = DP w + l=1 Vl ⊙ ⊤ Q (Q(Vl ⊙ w)) , where each matrix-vector product involving Q takes O((n + m)s) flops via the SVD of X, and thus computing Aw costs O(s2 (n+m)). Therefore, we can solve Aλ = b by the preconditioned conjugate gradient method, with Diag(diag(A)) as the preconditioner. Taking ξ = −∇f (X) for (A.2) yields the Gf (X) in Algorithm 1.
24
Y. YANG, B. GAO, AND Y.-X. YUAN REFERENCES
[1] P. Ablin and G. Peyré, Fast and accurate optimization on the orthogonal manifold without retraction, in International Conference on Artificial Intelligence and Statistics, PMLR, 2022, pp. 5636–5657. [2] P. Ablin, S. Vary, B. Gao, and P.-A. Absil, Infeasible deterministic, stochastic, and variance-reduction algorithms for optimization under orthogonality constraints, Journal of Machine Learning Research, 25 (2024), pp. 1–38. [3] P.-A. Absil, R. Mahony, and R. Sepulchre, Optimization Algorithms on Matrix Manifolds, Princeton University Press, 2008. [4] F. Andersson and M. Carlsson, Alternating projections on nontangential manifolds, Constructive Approximation, 38 (2013), pp. 489–525. [5] R. Andreani, E. G. Birgin, J. M. Martínez, and M. L. Schuverdt, On augmented Lagrangian methods with general lower-level constraints, SIAM Journal on Optimization, 18 (2008), pp. 1286–1309. [6] R. Andreani, K. R. Couto, O. P. Ferreira, and G. Haeser, Constraint qualifications and strong global convergence properties of an augmented Lagrangian method on Riemannian manifolds, SIAM Journal on Optimization, 34 (2024), pp. 1799–1825. [7] R. Andreani, K. R. Couto, O. P. Ferreira, G. Haeser, and L. F. Prudente, Global convergence of an augmented Lagrangian method for nonlinear programming via Riemannian optimization, SIAM Journal on Optimization, 36 (2026), pp. 466–501. [8] H. H. Bauschke, D. R. Luke, H. M. Phan, and X. Wang, Restricted normal cones and the method of alternating projections: applications, Set-Valued and Variational Analysis, 21 (2013), pp. 475–501. [9] H. H. Bauschke, D. R. Luke, H. M. Phan, and X. Wang, Restricted normal cones and the method of alternating projections: theory, Set-Valued and Variational Analysis, 21 (2013), pp. 431–473. [10] A. Beck and N. Hallak, On the minimization over sparse symmetric sets: projections, optimality conditions, and algorithms, Mathematics of Operations Research, 41 (2016), pp. 196–223. [11] R. Bergmann and R. Herzog, Intrinsic formulation of KKT conditions and constraint qualifications on smooth manifolds, SIAM Journal on Optimization, 29 (2019), pp. 2423– 2444. [12] N. Boumal, An Introduction to Optimization on Smooth Manifolds, Cambridge University Press, 2023. [13] N. Boumal, P.-A. Absil, and C. Cartis, Global rates of convergence for nonconvex optimization on manifolds, IMA Journal of Numerical Analysis, 39 (2019), pp. 1–33. [14] N. Boumal, B. Mishra, P.-A. Absil, and R. Sepulchre, Manopt, a Matlab toolbox for optimization on manifolds, The Journal of Machine Learning Research, 15 (2014), pp. 1455– 1459. [15] S. Budzinskiy, Quasioptimal alternating projections and their use in low-rank approximation of matrices and tensors, Numerische Mathematik, 157 (2025), pp. 1491–1535. [16] T. P. Cason, P.-A. Absil, and P. Van Dooren, Iterative methods for low rank approximation of graph similarity matrices, Linear Algebra and its Applications, 438 (2013), pp. 1863–1882. [17] S. Chen, Y. He, and W. Huang, Retractions by alternating projections, arXiv preprint arXiv:2605.17384, (2026). [18] S. Chen and W. Huang, Manifold identification and second-order algorithms for ℓ1 regularization on the Stiefel manifold. Talk at ICCOPT 2025, Los Angeles, CA. https: //iccopt2025usc.sched.com/event/21ZZR, 2025. [19] S. Chen, S. Ma, A. Man-Cho So, and T. Zhang, Proximal gradient method for nonsmooth optimization over the Stiefel manifold, SIAM Journal on Optimization, 30 (2020), pp. 210– 239. [20] M. Chu, N. Del Buono, L. Lopez, and T. Politi, On the low-rank approximation of data on the unit sphere, SIAM Journal on Matrix Analysis and Applications, 27 (2005), pp. 46–60. [21] D. Drusvyatskiy, A. D. Ioffe, and A. S. Lewis, Transversality and alternating projections for nonconvex sets, Foundations of Computational Mathematics, 15 (2015), pp. 1637–1651. [22] D. Drusvyatskiy and A. S. Lewis, Local linear convergence for inexact alternating projections on nonconvex sets, Vietnam Journal of Mathematics, 47 (2019), pp. 669–681. [23] H. Federer, Curvature measures, Transactions of the American Mathematical Society, 93 (1959), pp. 418–491.
OPTIMIZATION OVER THE INTERSECTION OF MANIFOLDS
25
[24] O. L. Frost, An algorithm for linearly constrained adaptive array processing, Proceedings of the IEEE, 60 (1972). [25] B. Gao, S. Vary, P. Ablin, and P.-A. Absil, Optimization flows landing on the Stiefel manifold, IFAC-PapersOnLine, 55 (2022), pp. 25–30. [26] F. Goyens and F. Feppon, The Riemannian landing method: from projected gradient flows to SQP, arXiv preprint arXiv:2603.24309, (2026). [27] L. Hörmander, The Analysis of Linear Partial Differential Operators III, Springer, Berlin, 1985. [28] A. D. Ioffe, Transversality in variational analysis, Journal of Optimization Theory and Applications, 174 (2017), pp. 343–366. [29] P. Jawanpuria, M. Meghwanshi, and B. Mishra, Low-rank approximations of hyperbolic embeddings, in 2019 IEEE 58th conference on decision and control, IEEE, 2019, pp. 7159– 7164. [30] X. Jia, C. Kanzow, P. Mehlitz, and G. Wachsmuth, An augmented Lagrangian method for optimization problems with structured geometric constraints, Mathematical Programming, 199 (2023), pp. 1365–1415. [31] Z. Lai and A. Yoshise, Riemannian interior point methods for constrained optimization on manifolds, Journal of Optimization Theory and Applications, 201 (2024), pp. 433–469. [32] J. M. Lee, Smooth Manifolds, Springer, 2012. [33] E. Levin, J. Kileel, and N. Boumal, The effect of smooth parametrizations on nonconvex optimization landscapes, Mathematical Programming, 209 (2025), pp. 63–111. [34] A. S. Lewis, D. R. Luke, and J. Malick, Local linear convergence for alternating and averaged nonconvex projections, Foundations of Computational Mathematics, 9 (2009), pp. 485–513. [35] A. S. Lewis and J. Malick, Alternating projections on manifolds, Mathematics of Operations Research, 33 (2008), pp. 216–234. [36] X. Li and Z. Luo, Normal cones intersection rule and optimality analysis for low-rank matrix optimization with affine manifolds, SIAM Journal on Optimization, 33 (2023), pp. 1333– 1360. [37] X. Li, N. Xiu, and S. Zhou, Matrix optimization over low-rank spectral sets: stationary points and local and global minimizers, Journal of Optimization Theory and Applications, 184 (2020), pp. 895–930. [38] C. Liu and N. Boumal, Simple algorithms for optimization on Riemannian manifolds with constraints, Applied Mathematics & Optimization, 82 (2020), pp. 949–981. [39] X. Liu, N. Xiao, and Y.-x. Yuan, A penalty-free infeasible approach for a class of nonsmooth optimization problems over the Stiefel manifold, Journal of Scientific Computing, 99 (2024), p. 30. [40] G. A. Miller, WordNet: a lexical database for English, Communications of the ACM, 38 (1995), pp. 39–41. [41] V. Neumann, Functional operators, The Geometry of Orthogonal Spaces, (1950). [42] M. Nickel and D. Kiela, Poincaré embeddings for learning hierarchical representations, Advances in neural information processing systems, 30 (2017). [43] J. Nocedal and M. L. Overton, Projected Hessian updating algorithms for nonlinearly constrained optimization, SIAM Journal on Numerical Analysis, 22 (1985), pp. 821–850. [44] J. Nocedal and S. J. Wright, Numerical Optimization, Springer, 2006. [45] D. Noll and A. Rondepierre, On local convergence of the method of alternating projections, Foundations of Computational Mathematics, 16 (2016), pp. 425–455. [46] M. Obara, T. Okuno, and A. Takeda, Sequential quadratic optimization for nonlinear optimization problems on Riemannian manifolds, SIAM Journal on Optimization, 32 (2022), pp. 822–853. [47] G. Olikier, P. Mlinarić, P.-A. Absil, and A. Uschmajew, The tangent cone to the real determinantal variety: various expressions and a proof, Set-Valued and Variational Analysis, 34 (2026), p. 8. [48] V. Ozolin, š, R. Lai, R. Caflisch, and S. Osher, Compressed modes for variational problems in mathematics and physics, Proceedings of the National Academy of Sciences, 110 (2013), pp. 18368–18373. [49] R. Peng, C. Zhu, B. Gao, X. Wang, and Y.-x. Yuan, Normalized tensor train decomposition, arXiv preprint arXiv:2511.04369, (2025). [50] R. T. Rockafellar, Augmented Lagrangians and applications of the proximal point algorithm in convex programming, Mathematics of operations research, 1 (1976), pp. 97–116. [51] R. T. Rockafellar and R. J.-B. Wets, Variational Analysis, vol. 317, Springer Science & Business Media, 2009.
26
Y. YANG, B. GAO, AND Y.-X. YUAN
[52] J. B. Rosen, The gradient projection method for nonlinear programming. Part I. Linear constraints, Journal of the Society for Industrial and Applied Mathematics, 8 (1960), pp. 181– 217. [53] J. B. Rosen, The gradient projection method for nonlinear programming. Part II. Nonlinear constraints, Journal of the Society for Industrial and Applied Mathematics, 9 (1961), pp. 514–532. [54] S. Schechtman, D. Tiapkin, M. Muehlebach, and E. Moulines, Orthogonal directions constrained gradient method: from non-linear equality constraints to Stiefel manifold, in The Thirty Sixth Annual Conference on Learning Theory, PMLR, 2023, pp. 1228–1258. [55] A. Schiela and J. Ortiz, An SQP method for equality constrained optimization on Hilbert manifolds, SIAM Journal on Optimization, 31 (2021), pp. 2255–2284. [56] W. Si and J. Malick, A unified landing framework for equality-constrained optimization, (2026), https://hal.science/hal-05487561. [57] Y. Song, P. Li, B. Gao, and K. Yuan, Distributed retraction-free and communicationefficient optimization on the Stiefel manifold, in International Conference on Machine Learning, 2025. [58] B. Vandereycken, Low-rank matrix completion by Riemannian optimization, SIAM Journal on Optimization, 23 (2013), pp. 1214–1236. [59] S. Vary, P. Ablin, B. Gao, and P.-A. Absil, Optimization without retraction on the random generalized Stiefel manifold, in Proceedings of the 41st International Conference on Machine Learning, vol. 235, PMLR, 2024, pp. 49226–49248. [60] N. Xiao, T. Tang, S. Wang, and K.-C. Toh, An exact penalty approach for equality constrained optimization over a convex set, arXiv preprint arXiv:2505.02495, (2025). [61] N. Xiao, S. Wang, T. Tang, and K.-C. Toh, A quadratically convergent alternating projection method for nonconvex sets, arXiv preprint arXiv:2511.22916, (2025). [62] X. Xiong, B. Gao, and P.-A. Absil, A second-order method landing on the Stiefel manifold via Newton–Schulz iteration, arXiv preprint arXiv:2605.02838, (2026). [63] W. H. Yang, L.-H. Zhang, and R. Song, Optimality conditions for the nonlinear programming problems on Riemannian manifolds, Pacific Journal of Optimization, 10 (2014), pp. 415–434. [64] Y. Yang, B. Gao, and Y.-x. Yuan, Variational analysis of determinantal varieties, arXiv preprint arXiv:2511.22613, (2025). [65] Y. Yang, B. Gao, and Y.-x. Yuan, A space-decoupling framework for optimization on bounded-rank matrices with orthogonally invariant constraints, Mathematical Programming, (2026), pp. 1–53. [66] Y.-x. Yuan, A null space algorithm for constrained optimization, Advances in Scientific Computing, Science Press, Beijing, (2001), pp. 210–218. [67] Y. Zhou, C. Bao, C. Ding, and J. Zhu, A semismooth Newton based augmented Lagrangian method for nonsmooth optimization on matrix manifolds, Mathematical Programming, 201 (2023), pp. 1–61.