# The Riemannian Landing Method: from projected gradient flows to SQP

Florentin Goyens, Florian Feppon

► To cite this version:

Florentin Goyens, Florian Feppon. The Riemannian Landing Method: from projected gradient flows to SQP. 2026. hal-05460164

HAL Id: hal-05460164

https://hal.science/hal-05460164v1

Preprint submitted on 15 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
THE RIEMANNIAN LANDING METHOD:
FROM PROJECTED GRADIENT FLOWS TO SQP

FLORENTIN GOYENS¹*, FLORIAN FEPPON²

¹ ICTEAM Institute, UCLouvain, Louvain-la-Neuve, Belgium

² NUMA Unit, Department of Computer Science, KU Leuven, Belgium.

**ABSTRACT.** Landing methods have recently emerged in Riemannian matrix optimization as efficient schemes for handling nonlinear equality constraints without resorting to costly retractions. These methods decompose the search direction into tangent and normal components, enabling asymptotic feasibility while maintaining inexpensive updates. In this work, we provide a unifying geometric framework that reveals the landing algorithm to encompass several classical optimization methods, under suitable choices of Riemannian metric, such as projected and null-space gradient flows, Sequential Quadratic Programming (SQP), and a certain form of the Augmented Lagrangian method. In particular, we show that a quadratically convergent landing method essentially reproduces the quadratically convergent SQP method. These connections also allow us to propose a globally convergent method using adaptive step sizes. The backtracking linesearch satisfies an Armijo condition on a merit function, and does not require a priori knowledge of Lipschitz constants.

Our second key contribution is to analyze landing methods through a geometric parameterization of the metric in terms of fields of oblique projectors and associated metric restrictions. This viewpoint disentangles the roles of orthogonality, tangent, and normal metrics, and elucidates how to design the metric so as to obtain explicit tangent and normal updates. For matrix optimization, this framework not only recovers recent constructions in the literature for problems with orthogonality constraints, but also provides systematic guidelines for designing new metrics that admit closed-form search directions.

**Keywords.** Landing algorithm, Constrained nonlinear optimization, Stiefel manifold, Projected gradient flows, Sequential Quadratic Programming, Oblique projections, Differential geometry in optimization.

**AMS 2020 Subject classifications.** 49M05, 65K10, 49M29, 90C30, 53C21

CONTENTS

<table><tr><td>1. Introduction</td><td>2</td></tr><tr><td>2. Historical developments of landing-related methods</td><td>3</td></tr><tr><td>3. Riemannian geometry and notation conventions</td><td>5</td></tr><tr><td>Metric pseudoinverse and formulas for g-orthogonal projectors</td><td>7</td></tr><tr><td>4. Alternative definitions of the landing algorithm and metric construction via normal bundles.</td><td>7</td></tr><tr><td>4.1. Construction of the ambient metric g from a given normal bundle</td><td>8</td></tr><tr><td>4.2. Equivalence between landing algorithms and projected gradient flows</td><td>9</td></tr><tr><td>5. Links between landing methods and existing constrained optimization algorithms</td><td>11</td></tr><tr><td>5.1. Interpretation of the landing algorithm as an augmented Lagrangian method</td><td>11</td></tr><tr><td>5.2. Sequential Quadratic Programming is a particular case of the landing algorithm</td><td>12</td></tr><tr><td>5.3. Link with the Riemannian Newton method: all roads lead to SQP</td><td>13</td></tr><tr><td>6. Global convergence using adaptive step sizes</td><td>17</td></tr><tr><td>6.1. Assumptions</td><td>17</td></tr><tr><td>6.2. Line search procedure</td><td>18</td></tr><tr><td>7. Geometric design of tangent and normal terms for matrix optimization with orthogonality constraints</td><td>22</td></tr><tr><td>7.1. Landing algorithm for the Euclidean metric</td><td>23</td></tr><tr><td>7.2. Metric design based on an explicit choice of projection operators</td><td>24</td></tr><tr><td>References</td><td>29</td></tr></table>

Date: January 15, 2026.

* Corresponding author. Email: goyensflorentin@gmail.com.

**Acknowledgements.** This work was supported by the Fonds de la Recherche Scientifique-FNRS under Grant no T.0001.23. F. Feppon was supported by the Flanders Research Foundation (FWO) under Grant G001824N.
# 1. INTRODUCTION

The optimization of functions subject to nonlinear constraints is a central problem across engineering, scientific computing, and machine learning. In recent years, the concept of landing methods has emerged in the Riemannian matrix optimization community after the works of Ablin and Peyré (2022) and Gao et al. (2022). The essential idea of the landing algorithm is to decompose the search direction

$$ (1.1) \qquad d = d_T + d_N $$

into a tangent component $d_T$, which is responsible for decreasing the objective function without worsening the constraint violation, and a normal component $d_N$, whose role is to progressively drive iterates toward the constraint manifold by reducing infeasibility.

In the context of matrix optimization, this approach enables one to enforce orthogonality constraints asymptotically, thus avoiding the computationally expensive retraction operations— i.e., nonlinear projections onto the feasible set—required by standard feasible Riemannian algorithms (Absil et al., 2008; Absil and Malick, 2012). Since evaluating the normal term $d_N$ is typically computationally cheaper than computing retractions, the method is appealing and makes geometric optimization potentially viable for large-scale problems, such as imposing the orthogonality of weight matrices during the training of Large Language Models (LLMs) or other resource-intensive machine learning tasks (Hu et al., 2025).

Consider an equality-constrained optimization problem of the form

$$ (P) \qquad \begin{cases} \min_{x \in \mathcal{E}} & f(x) \\ \text{s.t.} & c(x) = 0, \end{cases} $$

where $f: \mathcal{E} \to \mathbb{R}$ and $c: \mathcal{E} \to \mathcal{F}$ are continuously differentiable—possibly nonconvex—mappings between finite-dimensional vector spaces $\mathcal{E}$ and $\mathcal{F}$ with $\dim(\mathcal{F}) < \dim(\mathcal{E})$. In its most general formulation, the landing algorithm can be written as

$$ (1.2) \qquad x_{k+1} = x_k + \alpha_k (d_T(x_k) + d_N(x_k)), $$

where $\alpha_k > 0$ is an adaptive step size, and

$$ (1.3) \qquad d_T(x) = -\operatorname{grad}_{\mathcal{M}_x}^g f(x), $$

$$ (1.4) \qquad d_N(x) = -\nabla_g \psi(x) \text{ with } \psi(x) = \|c(x)\|_{\mathcal{F}}^2/2. $$

The tangent and normal directions $d_T(x)$ and $d_N(x)$ correspond respectively to the Riemannian constrained gradient of $f$ on the level curve

$$ (1.5) \qquad \mathcal{M}_x := \{y \in \mathcal{E} : c(y) = c(x)\}, $$

and the unconstrained Riemannian gradient of the infeasibility measure $\psi(x) := \|c(x)\|_{\mathcal{F}}^2/2$. The role of the tangent term is to decrease the objective function $f$ without increasing the violation of the constraints, while the purpose of the normal term is to ‘land’ the iterates $x_k$ back toward feasibility by descending along the constraint violation $c(x)$. Importantly, these gradients are computed with respect to a (user-defined) metric field $g$ that endows $\mathcal{E}$ with a Riemannian structure. When $\mathcal{E} = \mathbb{R}^n$ and $\mathcal{F} = \mathbb{R}^m$ are equipped with the standard Euclidean inner product, one obtains $d_N(x) = -\nabla c(x) c(x)$ and

$$ (1.6) \qquad d_T(x) = (I_n - \nabla c(x)^\top (\nabla c(x) \nabla c(x)^\top)^{-1} \nabla c(x)) \nabla f(x), $$

the orthogonal projection of $\nabla f(x)$ onto the tangent space of $\mathcal{M}_x$.

The first contribution of this paper is to show that, under suitable choices of the metric field $g$, the landing method (1.2) is equivalent to existing algorithms that have emerged under different names in other areas of optimization, such as the Corrective Gradient Projection method (Frost, 1972; Gangi and Byun, 1976), projected gradient flows (Tanabe, 1980; Yamashita, 1980; Schropp and Singer, 2000), and null-space gradient flows (Feppon et al., 2020a; Feppon, 2024). A more detailed account of the historical development of methods exploiting the decomposition of the search direction into tangent and normal components is given below. In fact, we show in section 4.2 that the scheme (1.2) includes, as particular cases, the Sequential Quadratic Programming (SQP) algorithm (Boggs and Tolle, 1995) in its basic form (without trust regions), as well as an Augmented Lagrangian method (Tapia, 1977) with a particular update rule for the multipliers. Although the introduction of landing methods for matrix optimization can thus be viewed as a resurgence of older classical methods, the questions of how to choose the metric $g$ and how to derive global convergence guarantees for the iterative scheme (1.2), rather than for its interpretation as the discretization of an Ordinary Differential Equation (ODE), are relatively recent and open new perspectives on existing algorithms.
Building on these connections, we propose an adaptive procedure for selecting the step size $\alpha_k$, based on an Armijo condition for the decrease of a merit function, that guarantees global convergence of the algorithm (1.2) toward a feasible point satisfying first-order optimality conditions. The procedure and its analysis rely heavily on standart techniques for the globalization of SQP (Nocedal and Wright, 2006; Curtis et al., 2024). Prior works proving the asymptotic convergence of the discrete iterates assumed constant step sizes, where global convergence is obtained for a step size smaller than some unknown constant (Schechtman et al., 2023; Zhang et al., 2024; Song et al., 2025). The present work thus fills an important gap, since landing schemes with constant step sizes are very sensitive to parameter tuning. Incidentally, a landing scheme with adaptive step sizes has been recently proposed in the concurrent work (Shi and Wang, 2025), but it relies on Lipschitz constants, while our linesearch does not rely on any problem constant.

Our second key contribution is to make explicit how (1.3) and (1.4) depend on the metric field $g$, and how to design it to (i) make (1.3) quadratically convergent, or (ii), to facilitate the computation of the tangent and normal steps (1.3) and (1.4).

Regarding the first goal, we realize that setting naturally the tangent term $d_T(x)$ to a Riemannian Newton step, the landing scheme (1.2) achieves quadratic convergence for unit step sizes only for a carefully designed normal step $d_N(x)$ such that (1.2) essentially reproduces SQP iterates (proposition 5.6). In reference to the leitmotiv “all roads lead to Newton” proclaimed in (Absil et al., 2009), we can thus state here that all roads lead to SQP for unfeasible constrained optimization methods.

Addressing the second goal leads us to introduce the following viewpoint to design the metric field $g$: first, associate a normal space $N_x M_x$ to every tangent space $T_x M_x$ of the constraint manifold $M_x$. Then, define the value of the restriction of the metric to the tangent and normal spaces (see (4.6)). Mathematically, this is requires specifying a family of oblique projectors $x \mapsto \text{Proj}_x$ on the tangent space at $x \in M_x$. The corresponding tangent and normal directions $d_T(x_k)$ and $d_N(x_k)$ can be expressed explicitly in terms of the oblique projector and these metric restrictions (proposition 4.3). In the context of Riemannian matrix optimization with orthogonality constraint, computing the tangent direction requires in general solving a linear system involving the inversion of $\nabla c(x)\nabla c(x)^\top$, which translates to solving Sylvester equations. Through a careful selection of the projector mapping and the tangent and normal metric restrictions, fully explicit formulas for $d_T(x)$ and $d_N(x)$ can be obtained, thereby enabling more efficient evaluations.

The article is organized as follows. Section 2 provides a detailed account of historical developments of constrained optimization methods closely related to the landing approach.

Section 3 reviews the background in differential geometry needed to introduce the notation conventions adopted in this work.

Section 4 establishes a formal link between the landing method and earlier projected gradient flows approaches. Such is achieved by exploiting the parameterization of the metric via a field of oblique projectors, together with the tangent and normal metric restrictions.

Section 5 establishes new connections between the landing methods and existing algorithms. We first show that (1.2) encompasses a version of the Augmented Lagrangian method with least-squares multiplier update, before proving that the Sequential Quadratic Programming (SQP) method is a particular instance of the Riemannian landing method. Then, we deepen known connections between SQP and Riemannian Newton method, formally establishing that a locally quadratically convergent landing method (1.2) is locally quadratically convergent essentially reproduces quadratically convergent SQP iterates, up to a negligible correction of the search direction.

In section 6, we establish global convergence of the Riemannian Landing Method in its most general form (1.2), with adaptive step sizes. The framework adapts well-known ideas for the gloablization of SQP methods based on a merit function.

Finally, section 7 presents several metric designs for optimization under orthogonality constraints, guided by a natural explicit choice of oblique projectors onto the tangent spaces of the constraint manifold. We detail the steps leading to explicit formulas for the normal and tangent components $d_T(x)$ and $d_N(x)$. This viewpoint clarifies why the metric fields considered in recent publications, e.g., (Ablin and Peyré, 2022; Goyens et al., 2026), lead to explicit formulas for $d_T(x)$ and $d_N(x)$, and it also allows us to determine how to design new metrics possessing these properties.
## 2. HISTORICAL DEVELOPMENTS OF LANDING-RELATED METHODS

To the best of our knowledge, the first trace of an algorithm decomposing the search direction as in (1.1) is the Gradient Projection Method formalized by Rosen (1960, 1961). For nonlinear constraints, this method projects at every step the gradient of the objective function onto the tangent constraint hyperplane—the null space of the linearized constraints—using (1.6), or onto the tangent cone in the case of inequality constraints. A similar algorithm was independently devised by Booker and Ong (1971) for geophysics signal processing with equality constraints. Shortly thereafter, the convergence rate of the method, considering geodesic retractions and equality constraints, was formally analyzed by Luenberger (1972).

However, finite steps on curved manifolds induce numerical infeasibility, which raised the need for an explicit correction term. This led to the development of the Corrective Gradient Projection (CGP) method (Frost, 1972; Gangi and Byun, 1976), which augments the tangent direction with a normal component correcting deviations from the constraints. This normal component,

$$ (2.1) \qquad d_N(x) = -\nabla c(x)^{\top}(\nabla c(x)\nabla c(x)^{\top})^{-1}c(x), $$

acts as a ‘pseudoinverse’ or ‘Newton-like’ step for finding a zero of the constraint.

In parallel with the development of CGP, Tanabe (1974) proposed to interpret Rosen’s gradient projection method as the discretization of the dynamical system $\dot{x} = d_T(x)$. Later, Yamashita (1980) introduced a corrected version of this projected gradient flow,

$$ (2.2) \qquad \begin{aligned} \dot{x} &= d_T(x) + d_N(x) \\ &= -(I_n - \nabla c(x)^{\top}(\nabla c(x)\nabla c(x)^{\top})^{-1}\nabla c(x))\nabla f(x) - \nabla c(x)^{\top}(\nabla c(x)\nabla c(x)^{\top})^{-1}c(x), \end{aligned} $$

rediscovering the pseudoinverse step of Frost (1972). This flow was later considered by Evtushenko and Zhadan (1994) to develop a barrier method for linear programming. Subsequently, Schropp and Singer (2000) proved that any bounded solution of (2.2) converges to a steady-state equilibrium, and stable equilibria coincide with local minimizers of (P). Schropp and Singer (2000) suggested that inequality constraints can be addressed with (2.2) through equalization via the introduction of additional slack variables; however, this produces an exponential growth in false equilibrium points as the number of inactive constraints increases (Jongen and Stein, 2003). To overcome this, Jongen and Stein (2004) and Shikhman and Stein (2009) proposed variations of (2.2) (without the corrective term), yielding modified gradient flows that implicitly incorporates the effects of slack variables.

A distinct differential-equation-based approach for problems with equality and inequality constraints emerged in shape and topology optimization under the name null space gradient flows or Null Space Optimizer (Feppon et al., 2020a). This method extends the projected flow (2.2) to general constrained problems and is compatible with the infinite-dimensional setting of shape optimization based on Hadamard’s boundary variation method (Henrot and Pierre, 2018; Allaire et al., 2021). This compatibility is possible thanks to the central role of the metric in converting differential into gradients and the interpretation of shape updates as retractions on an abstract manifold. Its tangent direction is the projection of the objective gradient onto the tangent cone of active or violated constraints, akin to the gradient projection of Rosen (1961), here computed by solving a quadratic program. The absence of hyperparameters makes the method particularly effective for topology and shape optimization (Feppon et al., 2020b). An alternative construction of tangent and normal terms for equality-constrained shape optimization was proposed in Wegert et al. (2023). A variant of the Null Space Optimizer for dealing with constraints with sparse Jacobian matrix was later introduced in Feppon (2024). Before these lines of work, let us mention that related gradient projection ideas had appeared in topology optimization: Yulin and Xiaoming (2004) projected shape gradients on tangent cones, while Barbarosie et al. (2020) proposed a method close to the Null Space Optimizer for equality constraints, together with a heuristic active-set strategy to handle inequalities.

More recently, the landing method was introduced in Ablin and Peyré (2022); Gao et al. (2022) for matrix optimization with orthogonality constraints without the need for retractions. This contrasts with the so far predominant Riemannian optimization methods, such as gradient flows or Riemannian Newton methods on the constraint manifold, which leverage intrinsic differential geometry and retractions Absil et al. (2008); Absil and Malick (2012).

Although the landing method introduced in Ablin and Peyré (2022) turns out to be very close in spirit with the projected gradient flow (2.2) of Yamashita (1980), there is a key difference: the normal term
in (1.2) is not a pseudoinverse step, but the unconstrained gradient of a penalty function. However, we show below that both coincide upon a suitable metric choice (proposition 4.4).

The landing algorithm was extended to general equality constraints and stochastic settings by Schechtman et al. (2023), who renamed it as *Orthogonal Directions Constrained Gradient Method*. The authors prove optimal complexity guarantees for the discretized algorithm, namely $O(\epsilon^{-2})$ iterations to reach an $\epsilon$-KKT point. To our knowledge, this is the first convergence proof for the discrete algorithm (1.2); earlier guarantees concerned the continuous system (2.2), assuming that the chosen discretization approximates the ODE trajectories sufficiently well. This analysis was further strengthened in Vary et al. (2024) in stochastic settings, using a smooth merit function to prove global convergence even in the presence of non-tangent random errors.

More recently, Zhang et al. (2024) applied the landing method to accelerate low-rank adaptation in fine-tuning large language models, proving convergence on the Stiefel manifold using only a sufficiently small constant step size. Goyens et al. (2026) derived landing flows on the Stiefel manifold for a family of metrics extending the $\beta$-metric of Hüper et al. (2021); Mataigne et al. (2025). Song et al. (2025) adapted the landing algorithm for distributed stochastic optimization on the Stiefel manifold. Lately, Shi and Wang (2025) proposed an extension called *Adaptive Directional Decomposition Methods* to incorporate inequality constraints. The tangent term is obtained by solving a quadratic program, similarly to the null-space gradient-flow strategy of Feppon et al. (2020a) for equality and inequality constrained problems. They show global convergence guarantees using an adaptive step size rule, which relies on the knowledge of Lipschitz constants.

## Links with null-space projections and SQP

The use of null-space projections to compute descent directions on the constraint manifold is certainly universal. It is therefore not surprising that formulas analogous to (1.6) and (2.1) have appeared in numerous alternative optimization frameworks, and even in robotics, where null-space control is used to exploit actuation redundancy and allow robots to perform multiple tasks simultaneously (Dietrich et al., 2015).

Most prominently, the decomposition (1.1) plays a central role in so-called *null-space methods*, where it is employed to isolate the objective descent within the null space and to construct reduced-Hessian approximations (Nocedal and Overton, 1985; Yuan, 2001; Nie, 2004; Biros and Ghattas, 2005; Nocedal and Wright, 2006; Berahas et al., 2024). These ideas lead to expressions directly comparable to the tangent and normal components in (2.2). Such decompositions into null-space and range-space steps have become standard in modern SQP implementations, most notably in the SNOPT software, where they are used both to reduce the dimension of the quadratic programming subproblem and to exploit reduced Hessians (Gill et al., 2002; Gill and Wong, 2012; Gill et al., 2015; Gharaei et al., 2023; Fang et al., 2024). They have also been exploited in interior-point methods (Liu and Yuan, 2010; Nocedal et al., 2014). Finally, we note that the least-squares multiplier $\lambda(x) = -(\nabla c(x)\nabla c(x)^\top)^{-1}\nabla c(x)\nabla f(x)$ is likewise employed in certain variants of the Augmented Lagrangian Method (Tapia, 1977; Conn et al., 1991), making these variants equivalent to landing-type algorithms (see proposition 5.2).

Finally, connections between SQP and differential-equation-based methods have also been pointed out repeatedly. For instance, Schropp and Singer (2000) showed that standard SQP schemes can be interpreted as variable step-size Euler discretizations of preconditioned versions of the flow (2.2). Precise relationships between the Riemannian Newton method and the quadratically convergent SQP method were established in Miller and Malick (2005) and Absil et al. (2009), and later exploited in Mishra and Sepulchre (2016) to design metric preconditioners accelerating the convergence of feasible Riemannian optimization methods on matrix manifolds. Recently, an analysis of the rate of convergence of the SQP algorithm with a constant quadratic term was derived by Bai and Mei (2018) through the interpretation as a Riemannian optimization method.

## 3. RIEMANNIAN GEOMETRY AND NOTATION CONVENTIONS

We consider the equality-constrained optimization problem (P) on the Euclidean spaces $\mathcal{E}$ and $\mathcal{F}$. We denote by $\langle \cdot, \cdot \rangle_{\mathcal{E}}$ and $\langle \cdot, \cdot \rangle_{\mathcal{F}}$ their respective inner products, and by $\|\cdot\|_{\mathcal{E}}$ and $\|\cdot\|_{\mathcal{F}}$ their associated norms. The feasible set is denoted by $\mathcal{M} := \{x \in \mathcal{E}: c(x) = 0\}$. In the applicable section 7 of this paper, $\mathcal{E}$ and $\mathcal{F}$ are thought of as matrix subspaces, but in the other parts of this section, the reader may think of these spaces as $\mathcal{E} \simeq \mathbb{R}^n$ and $\mathcal{F} \simeq \mathbb{R}^m$ through a proper vectorial indexing. If $V \subset \mathcal{E}$ is a subspace of $\mathcal{E}$, we denote by $V^{\perp, \mathcal{E}} \subset \mathcal{E}$ its orthogonal subspace with respect to the Euclidean metric of $\mathcal{E}$.
Throughout the paper, we denote by Dc the differential of the constraint function c: Dc(x) : E → F is the linear map such that

$$c(x + h) = c(x) + Dc(x)h + o(h) \text{ with } \frac{\|o(h)\|_{\mathcal{F}}}{\|h\|_{\mathcal{E}}} \to 0 \text{ as } \|h\|_{\mathcal{E}} \to 0.$$

We denote by D the set of points $x \in E$ at which the constraints satisfy the Linear Independence Constraint Qualification (LICQ) condition:

$$(3.1) \qquad \mathcal{D} = \{x \in E: \mathrm{rank}(Dc(x)) = m\}.$$

We recall that the set $M_x$ defined in (1.5) is a smooth manifold when it is included in D. Every set $M_x$ is a level curve of the constraint function c, called a layer manifold (Goyens et al., 2024). The tangent space of $M_x$ at $x \in D$ is the null space of the differential of the constraint:

$$(3.2) \qquad T_x M_x = \ker(Dc(x)).$$

## Riemannian metric and normal space

Let g be a Riemannian metric on E, that is a smooth family of inner products $x \mapsto g_x(\cdot, \cdot)$ on E. With a small abuse of notation, we drop the subscript x when referring to the metric. The metric is represented by a smoothly varying family of operators $G(x) : E \to E$, symmetric positive definite with respect to the Euclidean product $\langle \cdot, \cdot \rangle_E$, such that

$$(3.3) \qquad g(\xi, \zeta) = \langle G(x) \xi, \zeta \rangle_{\mathcal{E}}, \quad \text{for all } \xi, \zeta \in \mathcal{E}.$$

The corresponding norm is denoted by $\|u\|_g = \sqrt{g(u, u)}$ for any $u \in E$. We write occasionally $g^\mathcal{E} = \langle \cdot, \cdot \rangle_\mathcal{E}$ for the Euclidean metric.

The metric g allows to define several important objects for optimization. First, it determines the normal space $N_x^g$ to the manifold $M_x$ at every $x \in D$:

$$(3.4) \qquad N_x^g M_x := \{v \in E: g(v, \xi) = 0 \text{ for all } \xi \in T_x M_x\}.$$

Throughout the paper, we denote by Proj$_{x,g}$: E → T$_x$M$_x$ the orthogonal projection operator onto T$_x$M$_x \oplus N_x^g M$_x$. This operator is the unique linear operator satisfying

$$(3.5) \qquad g(\xi, v - \mathrm{Proj}_{x,g}(v)) = 0 \quad \text{for all } \xi \in T_x M_x \text{ and } v \in \mathcal{E},$$

see e.g., (Absil et al., 2008, (3.37)). In section 3, we will adopt a point of view where the tangent projector is defined before the metric g; in that case, we will denote it simply by Proj$_x$. The orthogonal projector associated to the Euclidean metric $g^\mathcal{E}$ will be occasionally denoted by $\Pi_x$.

The metric enables to define gradients and adjoints, whose definitions are recalled in the next paragraphs.

## Riemannian gradients on E and on M_x

Let Df(x) : E → R denote the differential of f at x ∈ D. Through the celebrated Riesz representation theorem, there are several ways to identify the linear form Df(x) to a vector playing the role of a gradient. First, for any x ∈ D, the unconstrained Euclidean gradient of f is the unique element denoted by $\nabla_{\mathcal{E}}f(x) \in \mathcal{E}$ such that

$$\langle \xi, \nabla_{\mathcal{E}} f(x) \rangle_{\mathcal{E}} = Df(x)[\xi] \quad \text{for all } \xi \in \mathcal{E}.$$

The unconstrained Riemannian gradient of f with respect to the metric g is written $\nabla_g f$ and is defined, for every $x \in D$, as the unique element $\nabla_g f(x) \in T_x D \simeq \mathcal{E}$ that satisfies

$$(3.6) \qquad g(\nabla_g f(x), \xi) = Df(x)[\xi] \quad \text{for all } \xi \in T_x D \simeq \mathcal{E}.$$

The Riemannian metric induces in turn a constrained Riemannian gradient on the manifold $M_x$, denoted by $\mathrm{grad}_{M_x}^g f$. Given $x \in D$, $\mathrm{grad}_{M_x}^g f(x)$ is the unique vector in $T_x M_x$ satisfying

$$(3.7) \qquad g(\mathrm{grad}_{M_x}^g f(x), \xi) = Df(x)[\xi], \quad \text{for all } \xi \in T_x M_x.$$

The unconstrained Euclidean and Riemannian gradients are related by the identity

$$(3.8) \qquad \nabla_g f(x) = G(x)^{-1} \nabla_{\mathcal{E}} f(x).$$

Moreover, the constrained Riemannian gradient is the orthogonal projection of the unconstrained Riemannian gradient on the tangent space $T_x M_x$:

$$(3.9) \qquad \mathrm{grad}_{M_x}^g f(x) = \mathrm{Proj}_{x,g}(\nabla_g f(x)).$$
## Metric adjoints

Given a linear operator $A : \mathcal{E} \to \mathcal{F}$, the adjoint with respect to the metric $g$ is defined as the unique linear operator $A^{*,g} : \mathcal{F} \to \mathcal{E}$ satisfying

$$ (3.10) \qquad \langle A\xi, y \rangle_{\mathcal{F}} = g(\xi, A^{*,g}y), \quad \text{for all } \xi \in \mathcal{E}, y \in \mathcal{F}. $$

We denote by $A^{*,\mathcal{E}}$ the adjoint operator with respect to the Euclidean metric, which satisfies:

$$ (3.11) \qquad \langle A\xi, y \rangle_{\mathcal{F}} = \langle \xi, A^{*,\mathcal{E}}y \rangle_{\mathcal{E}}, \quad \text{for all } \xi \in \mathcal{E}, y \in \mathcal{F}. $$

The Riemannian and the Euclidean adjoint operators are related through the identity

$$ (3.12) \qquad A^{*,g} = G(x)^{-1} A^{*,\mathcal{E}}. $$

In the particular case where $\mathcal{E} = \mathbb{R}^n$ and $\mathcal{F} = \mathbb{R}^m$, the differential $Dc(x)$ is represented by its Jacobian matrix

$$ (3.13) \qquad Dc(x) = \begin{pmatrix} \nabla c_1(x)^\top \\ \vdots \\ \nabla c_m(x)^\top \end{pmatrix} \in \mathbb{R}^{m \times n}, $$

and the Euclidean adjoint of any matrix $A \in \mathbb{R}^{m \times n}$ is given by the usual transpose $A^{*,\mathcal{E}} = A^\top$.

## Metric pseudoinverse and formulas for $g$-orthogonal projectors

An important operator in optimization algorithms is the mapping

$$ Dc(x) Dc(x)^{*,g} : \mathcal{F} \to \mathcal{F}, $$

which is self-adjoint with respect to the inner product on $\mathcal{F}$ and positive definite for $x \in \mathcal{D}$. The right $g$-pseudoinverse of $Dc(x)$ is defined as the operator

$$ (3.14) \qquad Dc(x)^{\dagger,g} := Dc(x)^{*,g} (Dc(x) Dc(x)^{*,g})^{-1} : \mathcal{F} \to \mathcal{E}, $$

which satisfies $Dc(x)Dc(x)^{\dagger,g} = \text{Id}_{\mathcal{F}}$. The $g$-orthogonal projectors onto $T_x M_x$ and $N_x M_x$ read then respectively

$$ (3.15) \qquad \text{Proj}_{x,g} = \text{Id}_{\mathcal{E}} - (Dc(x))^{\dagger,g} Dc(x), \quad \text{Proj}_{x,g}^{\perp} = (Dc(x))^{\dagger,g} Dc(x). $$

We note that the $g$-normal space is $N_x^g M_x = \text{Range}(Dc(x)^{*,g})$.

# 4. ALTERNATIVE DEFINITIONS OF THE LANDING ALGORITHM AND METRIC CONSTRUCTION VIA NORMAL BUNDLES.

This section establishes the equivalence between the landing algorithm (1.2) and one closer to the projected gradient flow (2.2). Before this, we observe that there is no gain of generality in considering two possibly different metrics for the tangent and normal terms.

**Lemma 4.1.** Consider the landing scheme (1.2) with the tangent and 'normal' terms calculated in two different metrics $g_1$ and $g_2$:

$$ (4.1) \qquad d_T(x) = -\text{grad}_{M_x}^{g_1} f(x), $$

$$ (4.2) \qquad d_N(x) = -\nabla_{g_2} \psi(x) \text{ with } \psi(x) = \|c(x)\|_{\mathcal{F}}^2/2. $$

There exists a common metric $g$ such that $d_T(x) = -\text{grad}_{M_x}^g f(x)$ and $d_N(x) = -\nabla_g \psi(x)$.

*Proof.* Consider the metric $g$ given by

$$ g(\xi, \zeta) := g_1(\text{Proj}_{x,g_2}\xi, \text{Proj}_{x,g_2}\zeta) + g_2(\text{Proj}_{x,g_2}^{\perp}\xi, \text{Proj}_{x,g_2}^{\perp}\zeta). $$

It is clear that for any $\xi, \zeta \in T_x M_x$, $g(\xi, \zeta) = g_1(\xi, \zeta)$, which implies $\text{grad}_{M_x}^{g_1} f(x) = \text{grad}_{M_x}^{g_2} f(x)$. On the other hand, $\nabla_{g_2} \Psi(x)$ is the unique vector in $\mathcal{E}$ satisfying $g_2(\nabla_{g_2} \Psi(x), \xi) = \langle Dc(x)\xi, c(x) \rangle_{\mathcal{F}}$ for all $\xi \in \mathcal{E}$. But for any $\xi \in \mathcal{E}$, $g(\nabla_{g_2} \Psi(x), \xi) = g_2(\nabla_{g_2} \Psi(x), \xi) = \langle Dc(x)\xi, c(x) \rangle_{\mathcal{F}}$. This implies $\nabla_{g_2} \Psi(x) = \nabla_g \Psi(x)$. $\square$
We note that the initial paper on the landing algorithm Gao et al. (2022) considered a Riemannian metric for the tangent term and the Euclidean metric for the computation of the normal term. The previous lemma shows that the metric in the normal space can always be redefined to infer both terms from a common metric.

The remainder of the section is structured as follows. We propose in section 4.1 a procedure for defining the ambient metric g by first specifying the orthogonality, in other words, the normal bundle of $M_x$ for every $x \in \mathcal{E}$. This enables to explicit the dependence of the tangent and normal terms $d_T(x)$ and $d_N(x)$ to the orthogonal projectors on $T_x M_x$ and $N_x M_x$ and the restriction of the metric g to these spaces.

Then, we show the equivalence between the landing algorithm (1.2) and a variation of it, where the normal step is replaced with a pseudoinverse or ‘Newton’-like step of the form of (2.1), as considered in the projected gradient flow (2.2). More precisely, we show that the pseudoinverse step $d_N(x)$ depends only on the choice of the normal space and not on the restriction of the normal metric. Then, we prove that the normal metric can always be redefined so that the pseudoinverse step becomes the Riemannian gradient of the penalty functional $\psi(x) = ||c(x)||^2/2$ as in (1.4).

## 4.1. Construction of the ambient metric g from a given normal bundle

In the definitions (1.3) and (1.4), the metric g is defined a priori, which enables to compute the tangent and the normal terms from the formulas

$$ (4.3) \qquad d_T(x) = -\mathrm{Proj}_{x,g}(G(x)^{-1}\nabla_{\mathcal{E}}f(x)) $$

$$ (4.4) \qquad d_N(x) = -G(x)^{-1}\mathrm{Dc}(x)^{*,\mathcal{E}}c(x). $$

The first formula follows from (3.8), (3.9) and (3.12), while the second can be inferred from the identities (3.6) and (3.10):

$$ g(\nabla_g \psi(x), \xi) = \mathrm{D}\psi(x)[\xi] = \langle \mathrm{Dc}(x)\xi, c(x) \rangle_{\mathcal{F}} = g(\mathrm{Dc}(x)^{*g}c(x), \xi) = g(G(x)^{-1}\mathrm{Dc}(x)^{*\mathcal{E}}c(x), \xi) \quad \forall \xi \in \mathcal{E}. $$

An issue with formulas (4.3) and (4.4), in the context of Riemannian matrix optimization, is that closed-form expressions for that the projection mapping $\mathrm{Proj}_{x,g}$ or the inverse metric $G(x)^{-1}$ are not necessarily available.

Alternatively, the metric g can be constructed by first considering a smooth mapping of (oblique) linear tangent projectors

$$ x \mapsto \mathrm{Proj}_x, $$

that is, for every $x \in \mathcal{D}$, $\mathrm{Proj}_x$ is a given linear mapping on $\mathcal{E}$ satisfying

$$ \mathrm{Proj}_x \mathrm{Proj}_x = \mathrm{Proj}_x, \quad \mathrm{Range}(\mathrm{Proj}_x) = T_x M_x \quad \text{and} \quad \mathrm{Proj}_x = \mathrm{Id} \text{ on } T_x M_x. $$

This viewpoint is equivalent to attaching to every $x \in M_x$ a ‘normal’ space characterized by

$$ N_x M_x := \ker(\mathrm{Proj}_x). $$

Then, consider two symmetric operator fields

$$ G_T(x) : \mathcal{E} \to \mathcal{E}, \qquad G_N(x) : \mathcal{E} \to \mathcal{E}, $$

required to be positive definite on respectively $T_x M_x$ and $N_x M_x$: there exist $g_T, g_N > 0$ such that

$$ (4.5) \qquad \begin{aligned} \langle \xi, G_T(x)\xi \rangle_{\mathcal{E}} &\ge g_T ||\xi||_{\mathcal{E}}^2, & \forall \xi \in T_x M_x, \\ \langle \xi, G_N(x)\xi \rangle_{\mathcal{E}} &\ge g_N ||\xi||_{\mathcal{E}}^2, & \forall \xi \in N_x M_x. \end{aligned} $$

The ambient metric g can finally be defined as

$$ (4.6) \qquad g(\xi, \zeta) = \langle \xi, G(x)\zeta \rangle = \langle \mathrm{Proj}_x \xi, G_T(x)\mathrm{Proj}_x \zeta \rangle_{\mathcal{E}} + \langle \mathrm{Proj}_x^{\perp} \xi, G_N(x)\mathrm{Proj}_x^{\perp} \zeta \rangle_{\mathcal{E}}, $$

denoting by $\mathrm{Proj}_x^{\perp} := \mathrm{Id}_{\mathcal{E}} - \mathrm{Proj}_x$ the oblique projector on $N_x M_x$. In other words, $g(\xi, \zeta) = \langle \xi, G(x)\zeta \rangle_{\mathcal{E}}$ where $G(x) : \mathcal{E} \to \mathcal{E}$ is the block-diagonal operator

$$ (4.7) \qquad G(x) = \mathrm{Proj}_x^{*,\mathcal{E}} G_T(x) \mathrm{Proj}_x + (\mathrm{Proj}_x^{\perp})^{*,\mathcal{E}} G_N(x) \mathrm{Proj}_x^{\perp}. $$

Then, it is clear that:

(i) $G_T$ and $G_N$ are representers of the restriction of the metric to the tangent space $T_x M_x$ and the normal space $N_x M_x$;

(ii) $T_x M_x$ and $N_x M_x$ are g-orthogonal subspaces;
(iii) $\mathrm{Proj}_x = \mathrm{Proj}_{x,g}$ and $\mathrm{Proj}_x^\perp = \mathrm{Proj}_{x,g}^\perp$ are the $g$-orthogonal projectors on the decomposition $\mathcal{E} = T_x\mathcal{M}_x \oplus N_x\mathcal{M}_x$.

We further note that $\mathrm{Proj}_x^{*,\mathcal{E}}$ and $(\mathrm{Proj}_x^\perp)^*, \mathcal{E}$ are the linear projectors associated to the decomposition $\mathcal{E} = N_x\mathcal{M}_x^\perp, \mathcal{E} \oplus T_x\mathcal{M}_x^\perp, \mathcal{E}$, and that the operators

$$ (4.8) \qquad \begin{aligned} \widetilde{G}_T &:= \mathrm{Proj}_x^{*,\mathcal{E}} G_T(x) \mathrm{Proj}_x & : T_x\mathcal{M}_x \to (N_x\mathcal{M}_x)^\perp, \mathcal{E}, \\ \widetilde{G}_N &:= (\mathrm{Proj}_x^\perp)^*, \mathcal{E} G_N(x) \mathrm{Proj}_x^\perp & : N_x\mathcal{M}_x \to (T_x\mathcal{M}_x)^\perp, \mathcal{E}, \end{aligned} $$

are invertible due to (4.5). Denoting by $\widetilde{G}_T(x)^{-1}$ and $\widetilde{G}_N(x)^{-1}$ their inverses, the inverse of $G(x)$ reads

$$ (4.9) \qquad G(x)^{-1} = \mathrm{Proj}_x \widetilde{G}_T(x)^{-1} \mathrm{Proj}_x^{*,\mathcal{E}} + \mathrm{Proj}_x^\perp \widetilde{G}_N(x)^{-1} (\mathrm{Proj}_x^\perp)^*, \mathcal{E}. $$

**Remark 4.1.** It is clear that only the restrictions of $G_T$ and $G_N$ to respectively the tangent space $T_x\mathcal{M}_x$ and the normal space $N_x\mathcal{M}_x$ matter in the definition of $G(x)$; the definitions of $G_T(x)$ and $G_N(x)$ as operators on the whole set $\mathcal{E}$ is motivated by the need to define a differentiable metric field through (4.7).

**Remark 4.2.** This approach, whereby the differential structure of a manifold is specified by a mapping of projectors, has been considered to analyze the geometric structure of certain algebraic matrix decompositions, see Feppon and Lermusiaux (2019).

We can now express the tangent and normal steps (1.3) and (1.4) in terms of $\mathrm{Proj}_x$, $\widetilde{G}_T(x)^{-1}$ and $\widetilde{G}_N(x)^{-1}$.

**Proposition 4.2.** The tangent and normal space steps of (1.3) and (1.4) can be rewritten in terms of the projectors, the tangent and the normal metric $G_T(x)$ and $G_N(x)$ as

$$ (4.10) \qquad d_T(x) = -\widetilde{G}_T(x)^{-1} \mathrm{Proj}_x^{*,\mathcal{E}} \nabla_{\mathcal{E}} f(x), $$

$$ (4.11) \qquad d_N(x) = -\widetilde{G}_N(x)^{-1} \mathrm{Dc}(x)^*, \mathcal{E} c(x). $$

*Proof.* By (3.8), (3.9) and (4.9),

$$ \begin{aligned} d_T(x) &= -\mathrm{grad}_{\mathcal{M}_x}^g f(x) = -\mathrm{Proj}_x(G(x)^{-1} \nabla_{\mathcal{E}} f(x)) \\ &= -\mathrm{Proj}_x(\mathrm{Proj}_x \widetilde{G}_T(x)^{-1} \mathrm{Proj}_x^{*,\mathcal{E}} + \mathrm{Proj}_x^\perp \widetilde{G}_N(x)^{-1} (\mathrm{Proj}_x^\perp)^*, \mathcal{E}) \nabla_{\mathcal{E}} f(x) \\ &= -\mathrm{Proj}_x \widetilde{G}_T(x)^{-1} \mathrm{Proj}_x^{*,\mathcal{E}} \nabla_{\mathcal{E}} f(x) = -\widetilde{G}_T(x)^{-1} \mathrm{Proj}_x^{*,\mathcal{E}} \nabla_{\mathcal{E}} f(x). \end{aligned} $$

This proves (4.10). For (4.11), we write

$$ \begin{aligned} d_N(x) &= -\mathrm{Dc}^{*,g}(x)c(x) = -G(x)^{-1}\mathrm{Dc}^{*,\mathcal{E}}(x)c(x) = -G(x)^{-1}(\mathrm{Proj}_x^\perp)^*, \mathcal{E} \mathrm{Dc}^{*,\mathcal{E}}(x)c(x) \\ &= -\widetilde{G}_N(x)^{-1}\mathrm{Dc}(x)^*, \mathcal{E} c(x). \end{aligned} $$

□

**Remark 4.3.** Although formula (4.10) suggests that $d_T(x)$ could depend on the choice of the normal space through the projector $\mathrm{Proj}_x^{*,\mathcal{E}}$, this is not the case. This is visible in the definition (3.7) of the constrained Riemannian gradient. The operator $\widetilde{G}_T(x)^{-1} \mathrm{Proj}_x^{*,\mathcal{E}}$ depends thus solely on the restriction of the metric $g$ to the tangent space $T_x\mathcal{M}_x$, but not on the particular choice of normal space. In particular, the formula

$$ \widetilde{G}_T(x)^{-1} \mathrm{Proj}_x^{*,\mathcal{E}} = (\mathrm{Proj}_x^{*,\mathcal{E}} G_T(x) \mathrm{Proj}_x)^{-1} \mathrm{Proj}_x^{*,\mathcal{E}} $$

does not depend on the choice of projection projector $\mathrm{Proj}_x$ so long as $\mathrm{Range}(\mathrm{Proj}_x) = T_x\mathcal{M}_x$. In fact, one may verify that

$$ \widetilde{G}_T(x)^{-1} \mathrm{Proj}_x^{*,\mathcal{E}} \xi = \sum_{i=1}^{\dim(T_x\mathcal{M}_x)} \langle \xi, u_i \rangle \mathcal{E} u_i, $$

where $(u_i)_{1 \le i \le \dim(T_x\mathcal{M}_x)}$ is any $g$-orthonormal basis of $T_x\mathcal{M}_x$. However, (4.10) and (4.11) show that closed-form expressions are available for $d_T(x)$ and $d_N(x)$ if this is the case for (i) the projection operator $\mathrm{Proj}_x$ and (ii) the inverse of the restriction operators $\widetilde{G}_T(x)$ and $\widetilde{G}_N(x)$.
## 4.2. Equivalence between landing algorithms and projected gradient flows

Given an arbitrary smooth symmetric positive-definite operator field

$$H(x): \mathcal{F} \to \mathcal{F},$$

we consider the iterative scheme

$$(4.12) \qquad x_{k+1} = x_k - \alpha_k(d_T(x_k) + d_N(x_k)),$$

where the tangent term $d_T(x) = -\mathrm{grad}_{\mathcal{M}_x}^g f(x) \in T_x \mathcal{M}_x$ is unchanged while the normal vector field $d_N(x)$ is defined for $x \in \mathcal{D}$ by

$$(4.13) \qquad d_N(x) := -\mathrm{Dc}(x)^{\dagger,g} H(x) c(x) \in \mathbb{N}_x^g \mathcal{M}_x.$$

We refer to (4.13) as a 'pseudoinverse' step. The iterative algorithm (4.12) may be interpreted as the discretization of the ordinary differential equation

$$(4.14) \qquad \begin{aligned} \dot{x} &= -d_T(x) - d_N(x) \\ &= (\mathrm{Id}_\varepsilon - \mathrm{Dc}(x)^{*,g} (\mathrm{Dc}(x) \mathrm{Dc}(x)^{*,g})^{-1} \mathrm{Dc}(x)) \nabla_g f(x) \\ &\quad - \mathrm{Dc}(x)^{*,g} (\mathrm{Dc}(x) \mathrm{Dc}(x)^{*,g})^{-1} H(x) c(x), \end{aligned}$$

When $H(x) = \mathrm{Id}$ and $g = g^\varepsilon$ is the Euclidean metric, (4.14) coincides with the 'Newton-like' or 'pseudoinverse' step of the projected gradient flow (2.2).

Formula (4.13) offers a somewhat more explicit control of how optimization trajectories land on the constraint manifold than (1.4) based on the unconstrained Riemannian gradient. It is a descent direction for the infeasibility measure $\psi(x) = \|c(x)\|_{\mathcal{F}}^2/2$ with a descent value depending explicitly on $H(x)$ and independent of $g$:

$$(4.15) \qquad \begin{aligned} \mathrm{D}\psi(x)[d_N(x)] &= \langle \mathrm{Dc}(x)[d_N(x)], c(x) \rangle_{\mathcal{F}} \\ &= - \langle \mathrm{Dc}(x)\mathrm{Dc}(x)^{*,g} (\mathrm{Dc}(x)\mathrm{Dc}(x)^{*,g})^{-1} H(x) c(x), c(x) \rangle_{\mathcal{F}} \\ &= - \langle H(x)c(x), c(x) \rangle_{\mathcal{F}} < 0. \end{aligned}$$

More precisely, the normal step $d_N(x)$ is the minimum $g$-norm solution of the undetermined system $\mathrm{Dc}(x)[d_N] = -H(x)c(x)$:

$$(4.16) \qquad d_N(x) = \arg \min_{d_N \in \varepsilon} \frac{1}{2} \|d_N\|_g^2 \quad \text{such that} \quad \mathrm{Dc}(x)[d_N] = -H(x)c(x).$$

This implies that, at the continuous level, the constraint values along the gradient flow trajectories $x(t)$ solving (4.14) satisfy

$$\frac{d}{dt}c(x(t)) = -H(x(t))c(x(t)),$$

which entails $c(x(t)) = c(x(0)) \exp(-\int_0^t H(x(s))ds)$. Hence, the eigenvalues and the eigenvectors of the symmetric operator $H(x)$ determine the instantaneous decay rate of the components of the constraint violation vector $c(x) \in \mathcal{F}$.

When $H(x) = \mathrm{Dc}(x)\mathrm{Dc}(x)^{*,g}$, we retrieve $d_N(x) = -\mathrm{Dc}(x)^{*,g}c(x) = -\nabla_g \psi(x)$, the negative unconstrained Riemannian gradient of the infeasibility measure $\psi(x) = \|c(x)\|_{\mathcal{F}}^2/2$. The following propositions show that the converse result is true: firstly, (4.13) does not depend on the choice of the normal metric $G_N(x)$; secondly, it is actually possible to redefine the normal metric $G_N(x)$ without changing the values of $d_T(x)$ and $d_N(x)$ so that $d_N(x)$ becomes the negative Riemannian unconstrained gradient of $\psi$ in the new metric.

**Proposition 4.3.** The 'pseudoinverse' normal space step (4.13) depends only on the choice of orthogonal normal space $N_x \mathcal{M}_x$ through the normal projector $\mathrm{Proj}_x^\perp$. More precisely, it holds

$$(4.17) \qquad d_N(x) = -\mathrm{Proj}_x^\perp \mathrm{Dc}(x)^{\dagger, \varepsilon} H(x) c(x).$$

*Proof.* Recall that $\mathrm{Proj}_x^\perp = \mathrm{Dc}(x)^{\dagger,g} \mathrm{Dc}(x)$ (eq. (3.15)). Multiplying to the right by $\mathrm{Dc}(x)^{\dagger, \varepsilon}$ and using $\mathrm{Dc}(x)\mathrm{Dc}(x)^{\dagger, \varepsilon} = \mathrm{Id}_{\mathcal{F}}$, we obtain the identity $\mathrm{Dc}(x)^{\dagger,g} = \mathrm{Proj}_x^\perp \mathrm{Dc}(x)^{\dagger, \varepsilon}$ relating the $g$- and Euclidean right-pseudoinverses. Formula (4.17) follows by substituting into (4.13). $\square$
**Proposition 4.4.** The pseudoinverse normal step (4.13) (or (4.17)) can be rewritten as

$$ (4.18) \qquad d_N(x) = -\nabla_g \psi(x), $$

with $\psi(x) = ||c(x)||_F^2/2$ for the metric field $g$ defined through (3.3) and (4.7) by setting

$$ G_N(x) := Dc(x)^*, \mathcal{E} H(x)^{-1} Dc(x), $$

which defines a symmetric positive-definite operator on $N_x M_x$.

*Proof.* The operator $G_N(x)$ is clearly symmetric for the Euclidean inner product. To prove positive definiteness of $G_N(x)$ on $N_x^\mathcal{E} M_x$, consider any nonzero $\xi \in N_x M_x$. We find

$$ \langle G_N(x)\xi, \xi\rangle_\mathcal{E} = \langle Dc(x)^*, \mathcal{E} H(x)^{-1} Dc(x)\xi, \xi\rangle_\mathcal{E} = \langle H(x)^{-1} Dc(x)\xi, Dc(x)\xi\rangle_\mathcal{F} > 0, $$

since $\xi \in N_x^\mathcal{E} M_x$, $\ker Dc(x) \cap N_x = \{0\}$ and $H(x)^{-1}$ is symmetric definite positive. Thus, $G_N(x) \succ 0$ on $N_x M_x$ and (4.6) (leaving $G_T(x)$ unchanged) does define a metric $g$ on $\mathcal{E}$.

Let us recall (equation (4.11)) that the unconstrained Riemannian gradient of $\psi(x)$ is given by

$$ (4.19) \qquad \nabla_g \psi(x) = \widetilde{G}_N(x)^{-1} Dc(x)^*, \mathcal{E} c(x). $$

Let us express $\nabla_g \psi(x)$ in pseudoinverse form. Due to $\ker(Dc(x)) = T_x M_x$ and $\text{Range}(Dc(x)^*, \mathcal{E}) = T_x M_x^\perp, \mathcal{E}$, we have

$$ (4.20) \qquad \begin{aligned} Dc(x) \operatorname{Proj}_x^\perp &= Dc(x), & \operatorname{Proj}_x^*, \mathcal{E} Dc(x)^*, \mathcal{E} &= 0 \\ Dc(x) \operatorname{Proj}_x &= 0, & (\operatorname{Proj}_x^\perp)^*, \mathcal{E} Dc(x)^*, \mathcal{E} &= Dc(x)^*, \mathcal{E}. \end{aligned} $$

These identities imply

$$ \widetilde{G}_N(x) = (\operatorname{Proj}_x^\perp)^*, \mathcal{E} G_N(x) \operatorname{Proj}_x^\perp = (\operatorname{Proj}_x^\perp)^*, \mathcal{E} Dc(x)^*, \mathcal{E} H(x)^{-1} Dc(x) \operatorname{Proj}_x^\perp = Dc(x)^*, \mathcal{E} H(x)^{-1} Dc(x). $$

Right multiplying by $\operatorname{Proj}_x^\perp Dc(x)^\dagger, \mathcal{E}$, this leads to

$$ \begin{aligned} \widetilde{G}_N(x) \operatorname{Proj}_x^\perp Dc(x)^\dagger, \mathcal{E} &= Dc(x)^*, \mathcal{E} H(x)^{-1} Dc(x) \operatorname{Proj}_x^\perp Dc(x)^\dagger, \mathcal{E} \\ &= Dc(x)^*, \mathcal{E} H(x)^{-1} Dc(x) Dc(x)^\dagger, \mathcal{E} = Dc(x)^*, \mathcal{E} H(x)^{-1}, \end{aligned} $$

invoking the first line of (4.20). Left multiplying by $\widetilde{G}_N(x)^{-1} : T_x M_x^\perp, \mathcal{E} \to N_x M_x$, which is allowed because $\text{Range}(Dc(x)^*, \mathcal{E}) = T_x M_x^\perp, \mathcal{E}$, and right multiplying by $H(x)$, we obtain

$$ (4.21) \qquad \widetilde{G}_N(x)^{-1} Dc(x)^*, \mathcal{E} = \operatorname{Proj}_x^\perp Dc(x)^\dagger, \mathcal{E} H(x). $$

Combining (4.19) and (4.21) yields (4.18). $\square$

In other words, defining the normal term as a negative unconstrained Riemannian gradient (1.4) or as the pseudo inverse step (4.13) leads to the same family of optimization algorithms.

## 5. LINKS BETWEEN LANDING METHODS AND EXISTING CONSTRAINED OPTIMIZATION ALGORITHMS

In this section, we establish explicit connections between the landing algorithm and several well-known constrained optimization algorithms from the literature. Specifically, in section 5.1, we demonstrate that the landing algorithm, as defined by (1.2), can be interpreted as an Augmented Lagrangian Method with the Lagrange multiplier assigned the least-squares multiplier value at every iteration.

In section 5.2, we prove that the most basic version of Sequential Quadratic Programming (SQP), without trust-regions, is a particular case of the landing algorithm.

In section 5.3, we deepen the connection—previously observed by multiple authors, e.g. Miller and Malick (2005); Absil et al. (2009); Mishra and Sepulchre (2016)—between the quadratically convergent version of SQP and the Riemannian-Newton method on a manifold. We show that when the tangent step in (1.2) is chosen to be the Riemannian-Newton direction, the landing method with a “Newton-like” normal update $d_N(x) = -Dc(x)^\dagger g c(x)$ achieves quadratic convergence only for a distinguished choice of normal space, which precisely coincides with the one mandated by the SQP framework.
## 5.1. Interpretation of the landing algorithm as an augmented Lagrangian method

Consider the iterative scheme

$$ (5.1) \qquad x_{k+1} = x_k - \alpha_k \nabla_g \mathcal{L}_{\beta_k}(x_k, \lambda_k), $$

with the augmented Lagrangian $\mathcal{L}_{\beta_k}$ given by

$$ \mathcal{L}_{\beta_k}(x, \lambda) = f(x) + \langle \lambda, c(x) \rangle_{\mathcal{F}} + \frac{\beta_k}{2} \|c(x)\|_{\mathcal{F}}^2. $$

Here, $\beta_k$ is a penalty parameter which is not necessarily constant. In (5.1), the unconstrained Riemannian gradient $\nabla_g \mathcal{L}_{\beta_k}$ is taken with respect to the variable $x$.

**Proposition 5.1.** The Augmented Lagrangian iteration (5.1) coincides with the Riemannian landing iteration (4.12) with pseudoinverse step (4.13), setting $\lambda_k$ as the least-squares multiplier

$$ (5.2) \qquad \lambda_k := \arg \min_{\lambda \in \mathcal{F}} \| \nabla_g f(x_k) - Dc(x_k)^{*,g} \lambda \|_{\mathcal{F}}^2 = (Dc(x_k)Dc(x_k)^{*,g})^{-1} Dc(x_k) \nabla_g f(x_k). $$

$$ \text{and } H(x_k) := \beta_k Dc(x) Dc(x)^{*,g}. $$

*Proof.* The unconstrained Riemannian gradient of $\mathcal{L}_{\beta_k}$ with respect to $x$ at $x = x_k$ reads

$$ \begin{aligned} \nabla_g \mathcal{L}_{\beta_k}(x_k, \lambda_k) &= \nabla_g f(x_k) + Dc(x_k)^{*,g} \lambda_k + \beta_k Dc(x_k)^{*,g} c(x_k) \\ &= (\mathrm{Id}_{\mathcal{E}} - Dc(x_k)^{*,g} (Dc(x_k)Dc(x_k)^{*,g})^{-1} Dc(x_k)) \nabla_g f(x_k) + \beta_k Dc(x_k)^{*,g} c(x_k). \end{aligned} $$

Using (3.9) and (4.13) with $H(x) = \beta_k Dc(x) Dc(x)^{*,g}$, we obtain $\nabla_g \mathcal{L}_{\beta_k}(x_k, \lambda_k) = \mathrm{grad}_{\mathcal{M}_x}^g f(x_k) + d_N(x_k)$. $\square$

Thus, the landing algorithm is an augmented Lagrangian method with the Lagrange multipliers $\lambda_k$ being the least-squares multipliers, and for which the minimization of the augmented Lagrangian sub-problem consists in a single gradient step.

## 5.2. Sequential Quadratic Programming is a particular case of the landing algorithm

In its most basic form as described in Boggs and Tolle (1995), the Sequential Quadratic Programming (SQP) method considers the iterative sequence

$$ (5.3) \qquad x_{k+1} = x_k + \alpha_k d_k, $$

where $d_k$ is obtained by solving at each iteration the quadratic program

$$ (5.4) \qquad \begin{aligned} d_k := & \arg \min_{d \in \mathcal{E}} & f(x_k) + \langle d, \nabla_{\mathcal{E}} f(x_k) \rangle_{\mathcal{E}} + \frac{1}{2} \langle d, B_k d \rangle_{\mathcal{E}} \\ & \text{subject to} & Dc(x_k) d + c(x_k) = 0, \end{aligned} $$

given a sequence of symmetric positive definite operators $B_k : \mathcal{E} \to \mathcal{E}$. This method is locally convergent under mild assumptions on the selection of symmetric positive-definite matrices $B_k$. Setting $B_k = \nabla_{\mathcal{E}}^2 \mathcal{L}(x, \lambda_k)$ ($\nabla_{\mathcal{E}}^2$ being the Euclidean Hessian in the $x$ variable), with the Lagrangian

$$ \mathcal{L}(x, \lambda) := f(x) - \langle \lambda, c(x) \rangle_{\mathcal{F}}, $$

SQP achieves quadratic convergence around a KKT point for a unit step size when $\lambda_k \in \mathcal{F}$ is set to the Lagrange multiplier associated to the equality constraint of (5.4) of the previous iteration (see Nocedal and Wright (2006)). Globalization procedures also exist, based on merit functions or filters (Boggs and Tolle, 1995; Fletcher et al., 2002; Wächter and Biegler, 2005; Obara et al., 2022).

The following proposition shows that the SQP algorithm (5.3) forms a particular instance of the Riemannian landing method (1.2).

**Proposition 5.2.** Suppose that $B_k$ is a symmetric positive definite operator on $T_{x_k} \mathcal{M}_{x_k}$. Consider the space $\mathcal{N}_{x_k} \mathcal{M}_{x_k}$ defined as the $\mathcal{E}$-orthogonal space of $B_k T_{x_k} \mathcal{M}_{x_k}$:

$$ (5.5) \qquad \mathcal{N}_{x_k} \mathcal{M}_{x_k} := (B_k T_{x_k} \mathcal{M}_{x_k})^{\perp, \mathcal{E}}. $$

Let $g$ be the metric defined by

$$ g(\xi, \zeta) := \langle \mathrm{Proj}_{x_k} \xi, B_k \mathrm{Proj}_{x_k} \zeta \rangle_{\mathcal{E}} + \langle \mathrm{Proj}_{x_k}^{\perp} \xi, G_N(x_k) \mathrm{Proj}_{x_k}^{\perp} \zeta \rangle_{\mathcal{E}}, $$
where Proj$_{x_k}$ is the oblique projection on the decomposition $\mathcal{E} = \text{T}_{x_k} \mathcal{M}_{x_k} \oplus \text{N}_{x_k} \mathcal{M}_{x_k}$ and $\text{G}_N(x_k)$ any positive symmetric definite operator on $\text{N}_{x_k} \mathcal{M}_{x_k}$. Then, if $x_k \in \mathcal{D}$, the SQP direction $d_k$ of (5.4) is given by

$$ (5.6) \qquad d_k = -\operatorname{grad}_{\mathcal{M}_{x_k}}^g f(x_k) - \operatorname{Dc}(x_k)^{\dagger,g} c(x_k). $$

In other words, the SQP direction $d_k$ and the Riemannian landing direction $d_T(x_k) + d_N(x_k)$ based on the pseudoinverse formula with the choice $H(x_k) = \text{Id}_\mathcal{F}$ for the normal component (eq. (1.3) and (4.13)) coincide at every step k for the metric induced by $B_k$ on the tangent space $\text{T}_{x_k} \mathcal{M}_{x_k}$ and the normal space (5.5).

**Proof.** The space decomposition $\mathcal{E} = \text{T}_{x_k} \mathcal{M}_{x_k} \oplus \text{N}_{x_k} \mathcal{M}_{x_k}$ holds because $B_k$ is symmetric definite positive. Observe that any $d$ satisfying $\operatorname{Dc}(x_k)d + c(x_k)$ can be written as

$$ d = u + d_N(x_k), $$

with $d_N(x_k) = -\operatorname{Dc}(x_k)^{\dagger,g} c(x_k)$ and $u \in \text{T}_{x_k} \mathcal{M}_{x_k}$. The minimizer of (5.4) is $d_k = u_k + d_N(x_k)$, where $u_k$ is the solution to the quadratic unconstrained minimization problem

$$ (5.7) \qquad u_k = \arg \min_{u \in \text{T}_{x_k} \mathcal{M}_{x_k}} f(x_k) + \langle u + d_N(x_k), \nabla_{\mathcal{E}} f(x_k) \rangle_{\mathcal{E}} + \frac{1}{2} \langle u + d_N(x_k), B_k(u + d_N(x_k)) \rangle_{\mathcal{E}}. $$

Note further that for $u \in \text{T}_{x_k} \mathcal{M}_{x_k}$,

$$ \langle u, \nabla_{\mathcal{E}} f(x_k) \rangle_{\mathcal{E}} = \langle \text{Proj}_{x_k} u, \nabla_{\mathcal{E}} f(x_k) \rangle_{\mathcal{E}} = \langle u, \text{Proj}_{x_k}^{*,\mathcal{E}} \nabla_{\mathcal{E}} f(x_k) \rangle_{\mathcal{E}}, $$

and since $d_N(x_k) \in \text{N}_{x_k} \mathcal{M}_{x_k} = (B_k \text{T}_{x_k} \mathcal{M}_{x_k})^{\perp, \mathcal{E}}$,

$$ \begin{aligned} \langle u + d_N(x_k), B_k(u + d_N(x_k)) \rangle_{\mathcal{E}} &= \langle u, B_k u \rangle_{\mathcal{E}} + 2 \langle u, B_k d_N(x_k) \rangle_{\mathcal{E}} + \langle d_N(x_k), B_k d_N(x_k) \rangle_{\mathcal{E}} \\ &= \langle u, G_T(x_k) u \rangle_{\mathcal{E}} + 2 \langle d_N(x_k), B_k u \rangle_{\mathcal{E}} + \langle d_N(x_k), B_k d_N(x_k) \rangle_{\mathcal{E}} \\ &= \langle u, G_T(x_k) u \rangle_{\mathcal{E}} + \langle d_N(x_k), B_k d_N(x_k) \rangle_{\mathcal{E}}, \end{aligned} $$

where we denote $G(x_k) := B_k$. Consequently, by eliminating constant terms, the solution $u_k$ to (5.7) is also the minimizer of the quadratic problem

$$ (5.8) \qquad u_k = \arg \min_{u \in \text{T}_{x_k} \mathcal{M}_{x_k}} \langle u, \text{Proj}_{x_k}^{*,\mathcal{E}} \nabla_{\mathcal{E}} f(x_k) \rangle_{\mathcal{E}} + \frac{1}{2} \langle u, \widetilde{G}_T(x_k) u \rangle_{\mathcal{E}}, $$

where $\widetilde{G}_T(x_k) : \text{T}_{x_k} \mathcal{M}_{x_k} \to \text{N}_{x_k}^{\perp, \mathcal{E}}$ is the operator defined in (4.8). The solution to this quadratic program is $u_k = -\widetilde{G}_T(x_k)^{-1} \text{Proj}_{x_k}^{*,\mathcal{E}} \nabla_{\mathcal{E}} f(x_k)$, which is exactly $d_T(x_k)$ according to (4.10). $\square$

**Remark 5.1.** An alternative, more algebraic proof can be obtained if one assumes $B_k$ to be symmetric positive definite on the whole space $\mathcal{E}$ rather than only on $\text{T}_{x_k} \mathcal{M}_{x_k}$. In this case, the result is well-known and can be found stated in different forms in the optimization literature; see for instance (Boggs and Tolle, 1995, eq. (3.8)-(3.10)) or more recently (Feppon, 2024, Proposition 1). The idea of the proof is recalled here to highlight the link with the landing method (4.12). The KKT condition for (5.4) states that there exists a Lagrange multiplier $\lambda_k \in \mathcal{F}$ such that the minimizer $d_k$ reads

$$ (5.9) \qquad d_k = -B_k^{-1} \nabla_{\mathcal{E}} f(x_k) + B_k^{-1} \operatorname{Dc}(x_k)^{*,*} \lambda_k. $$

Inserting this expression into $\operatorname{Dc}(x_k)d_k = -c(x_k)$, we infer that

$$ \operatorname{Dc}(x_k) d_k = -c(x_k) = -\operatorname{Dc}(x_k) B_k^{-1} \nabla_{\mathcal{E}} f(x_k) + \operatorname{Dc}(x_k) B_k^{-1} \operatorname{Dc}(x_k)^{*,*} \lambda_k. $$

Since $\operatorname{Dc}(x_k)$ is full rank and $B_k^{-1}$ is assumed to be invertible, the operator $\operatorname{Dc}(x_k) B_k^{-1} \operatorname{Dc}(x_k)^{*,*}$ is symmetric positive definite and $\lambda_k$ is given by

$$ \lambda_k = -(\operatorname{Dc}(x_k) B_k^{-1} \operatorname{Dc}(x_k)^{*,*})^{-1} c(x_k) + (\operatorname{Dc}(x_k) B_k^{-1} \operatorname{Dc}(x_k)^{*,*})^{-1} \operatorname{Dc}(x_k) B_k^{-1} \nabla_{\mathcal{E}} f(x_k). $$

Substituting this expression into (5.9) yields

$$ (5.10) \qquad d_k = - \left( \operatorname{I}_n - B_k^{-1} \operatorname{Dc}(x_k)^{*,*} (\operatorname{Dc}(x_k) B_k^{-1} \operatorname{Dc}(x_k)^{*,*})^{-1} \operatorname{Dc}(x_k) \right) B_k^{-1} \nabla_{\mathcal{E}} f(x_k) \\ - B_k^{-1} \operatorname{Dc}(x_k)^{*,*} (\operatorname{Dc}(x_k) B_k^{-1} \operatorname{Dc}(x_k)^{*,*})^{-1} c(x_k), $$

which is (5.6) with the metric $g(\xi, \zeta) := \langle \xi, B_k \zeta \rangle$.
Conversely, given a metric g at $x_k$, the proof of proposition 5.2 shows that the “landing” step (5.6) can be obtained as the SQP step $d_k$ solving (5.4) by choosing

$$B_k := \operatorname{Proj}_{x_k}^{*,\mathcal{E}} G_T(x_k) \operatorname{Proj}_{x_k}.$$

Hence, the SQP algorithm (5.3) and (5.4) constitutes a strict subclass of landing algorithms (1.2) to (1.4). In SQP, the iterate is completely determined by the choice of the normal space and the tangent metric, whereas the landing framework additionally allows one to tune the normal metric for the normal step (1.4)—or, equivalently, to choose the operator H in (4.13).

Lastly, it is interesting to notice that, given a matrix $B_k$ generating the SQP direction $d_k$ through (5.4), the normal component can be modified without changing the tangent component by replacing $B_k$ with $\operatorname{Proj}_{x_k}^{*,\mathcal{E}} B_k \operatorname{Proj}_{x_k}$, where $\operatorname{Proj}_{x_k}$ is any linear projector on $T_{x_k}M_{x_k}$. After this change, the SQP direction reads (5.6) where g is a metric making $T_{x_k}M_{x_k}$ and $N_{x_k}M_{x_k} := \ker(\operatorname{Proj}_{x_k})$ g-orthogonal and satisfying $g(\xi, \zeta) = \langle \xi, B_k\zeta \rangle_{\mathcal{E}}$ for any $\xi, \zeta \in T_{x_k}M_{x_k}$.

## 5.3. Link with the Riemannian Newton method: all roads lead to SQP

We have recalled earlier that the SQP method (5.4) is quadratically convergent with $B_k = \nabla_{\mathcal{E}}^2 \mathcal{L}(x_k, \lambda_k)$ and $\lambda_k$ being the Lagrange multiplier of the quadratic program at the previous iteration. In fact, these settings make SQP equivalent to the Newton method for the functional $F(x, \lambda) := (\nabla_{\mathcal{E}} \mathcal{L}(x, \lambda), c(x))$ (Nocedal and Wright, 2006, Chapter 18). In view of proposition 5.2, this implies that in order to make the landing method

$$(5.11) \qquad x_{k+1} = x_k - \operatorname{grad}_{M_{x_k}}^g f(x_k) - Dc(x_k)^{\dagger,g} c(x_k)$$

locally quadratically convergent, it is sufficient to choose a metric g such that

(i) the normal space at $x_k$ is $N_{x_k}M_{x_k} := (\nabla_{\mathcal{E}}^2 \mathcal{L}(x_k, \lambda_k)T_{x_k}M_{x_k})^{\perp,\mathcal{E}}$;

(ii) the tangent metric is $G_T(x_k) = \nabla_{\mathcal{E}}^2 \mathcal{L}(x_k, \lambda_k)$.

On the other hand, it is known since Edelman et al. (1998); Absil et al. (2009) that when setting $G_T(x) = \nabla_{\mathcal{E}}^2 \mathcal{L}(x, \lambda^*(x))$ and $\operatorname{Proj}_x = \Pi_x$ the $\mathcal{E}$-orthogonal projection, the tangent term becomes a Riemannian Newton step:

$$-\operatorname{grad}_{M_x}^g f(x) = -\operatorname{Hess}_{M_x}^{\mathcal{E}} f(x)^{-1} \operatorname{grad}_{M_x}^{g_{\mathcal{E}}} f(x).$$

Here, $\operatorname{Hess}_{M_x}^{\mathcal{E}} f(x) = \Pi_x \nabla_{\mathcal{E}}^2 \mathcal{L}(x, \lambda(x)) \Pi_x$ is the Riemannian Hessian with respect to the Riemannian connection associated to the Euclidean metric; this follows from the fact that $\xi := \operatorname{grad}_{M_x}^g f(x)$ satisfies

$$\tilde{G}_T(x)\xi = \Pi_x \nabla_{\mathcal{E}}^2 \mathcal{L}(x, \lambda(x)) \Pi_x \xi = -\Pi_x \nabla_{\mathcal{E}} f(x) = -\operatorname{grad}_{M_x}^{g_{\mathcal{E}}} f(x).$$

In Absil et al. (2009), this property enabled to interpret the Feasibly Projected SQP (FP-SQP) method, i.e.

$$(5.12) \qquad x_{k+1} = R_{x_k}(\alpha_k d_k),$$

with $d_k$ being the SQP direction (5.4) and $B_k = \nabla_{\mathcal{E}}^2 \mathcal{L}(x_k, \lambda_k^*(x_k))$, as a Riemannian Newton method on the manifold $M$. Here, the multiplier is not the SQP multiplier $\lambda_k$ as above but the least-squares multiplier $\lambda^*(x_k)$ where

$$\lambda^*(x) = (Dc(x)Dc(x)^{*\mathcal{E}})^{-1}Dc(x)\nabla_{\mathcal{E}} f(x),$$

however, both choices are asymptotically equivalent since $\lambda_k \to \lambda^*(x^*)$ as $k \to +\infty$.

The FP-SQP method (5.12) combines the SQP step $d_k$ with $R_{x_k}$ a retraction on $M$ to ensure the feasibility of the iterates at every iteration.

In what follows, we deepen this connection, by showing that the conditions (i) and (ii) are in some sense necessary to make (5.11) locally quadratically convergent. Consider the iterative scheme

$$(5.13) \qquad x_{k+1} = x_k - \operatorname{Hess}_{M_{x_k}}^{g,\mathcal{A}} f(x_k)^{-1} \operatorname{grad}_{M_{x_k}}^g f(x_k) - Dc(x_k)^{\dagger,g} c(x_k),$$

which combines a Riemannian-Newton tangent step $-\operatorname{Hess}_{M_x}^{g,\mathcal{A}} f(x_k)^{-1} \operatorname{grad}_{M_{x_k}}^g f(x_k)$ on the manifold $M_{x_k}$ with the ‘Newton-like’ pseudoinverse normal step $-Dc(x_k)^{\dagger,g} c(x_k)$ aiming to reduce the constraint violation. The scheme (5.13) can be seen as a particular case of the landing method (4.12), up to redefining the metric in the tangent space.
Let us recall that the Riemannian Hessian is defined with respect to the metric $g$ and the choice of an affine connection $\nabla^A$ by the formula

$$ \mathrm{Hess}_{M_x}^{g,A} f(x) \xi := \nabla_{\xi}^A \mathrm{grad}_{M_x}^g f(x), $$

where an affine connection $\nabla^A$ is a bilinear operator on the tangent bundle that verifies Leibniz rule $(\nabla_{\xi}^A \zeta) \in T_x M_x$ and $\nabla_{\xi}^A (f\zeta) = D_{\xi}f\zeta + f\nabla_{\xi}^A \zeta$ for any smooth real function $f$ and any $\xi, \zeta \in T_x M_x$), see Absil et al. (2009), chapter 5. The superscript $\cdot^A$ in the notation $\nabla^A$ of the affine connection is suggested by the next lemma, which states that affine connections on a Riemannian submanifold can be parameterized by the choice of a tangent bilinear term $A(\xi, \zeta)$.

**Lemma 5.3.** *Any affine connection on $M_x$ can be written as*

$$ (5.14) \qquad \nabla_{\xi} \eta \equiv \nabla_{\xi}^{A} \eta := \mathrm{Proj}_{x}(D_{\xi} \eta) + A(\xi, \eta), $$

where $A(\xi, \eta) : T_x M_x \times T_x M_x \to T_x M_x$ is a vectorial bilinear form with values on the tangent space $T_x M_x$, and $\mathrm{Proj}_x$ is any projection operator on $T_x M_x$. Reciprocally, (5.14) defines an affine connection.

*Proof.* According to Gauss formula, any affine connection can be written as

$$ (5.15) \qquad \nabla_{\xi} \eta = d_{\xi} \eta - \Gamma(\xi, \eta), $$

where $\Gamma(\xi, \eta)$ is the Christoffel symbol of the connection. Since $\nabla_{\xi} \eta \in T_x M_x$, we must have

$$ \mathrm{Proj}_{x}^{\perp}(\Gamma(\xi, \eta)) = \mathrm{Proj}_{x}^{\perp}(D_{\xi}\eta), $$

so that $\Gamma(\xi, \eta) = \mathrm{Proj}_{x}^{\perp}(D_{\xi}\eta) - A(\xi, \eta)$ with $A(\xi, \eta) := -\mathrm{Proj}_{x}(\Gamma(\xi, \eta))$. Substituting this expression back into (5.15) yields (5.14). Conversely, we verify that for a given vector bilinear form $A(\xi, \eta)$, formula (5.14) satisfies all the axioms of an affine connection. $\square$

Recall that on the one hand, the scheme $x_{k+1} = x_k - Dc(x_k)^{\dagger,g}c(x_k)$ generates iterates that approach the feasible set at a quadratic rate, and on the other hand, the Riemannian-Newton method (5.12) with $d_k = -\mathrm{Hess}_{M_{x_k}}^{g,A} \mathrm{grad}_{M_{x_k}}^g$ is locally quadratically convergent to a critical point; these properties hold for any choice of metric $g$ and any choice of affine connection $\nabla^A$ (Absil et al., 2008). It could therefore be tempting to think that the scheme (5.13) combining this two directions should be automatically locally quadratically convergent near critical points of $f$ on the manifold $M$.

We show in the following that this is not true: the local quadratic convergence is obtained only by selecting a metric $g$ such that the normal space at $x_k$ is $N_{x_k} M_{x_k} = (\nabla_{\xi}^2 \mathcal{L}(x_k, \lambda^*(x_k))T_x M_x)^{\perp, \mathcal{E}}$ (or a convergent perturbation of it). This means that quadratic convergent iterates “land” on the manifold $M$ only through a preferred distinguished direction of the space, which is exactly the one taken by SQP.

We prove this result in two steps. First, we observe that the choice of affine connection influences the Hessian only up to a term vanishing near critical points.

**Proposition 5.4.** *There exists a bilinear tangent mapping $\tilde{A}_g : T_x M_x \times T_x M_x \to T_x M_x$ continuous with respect to $g$ and $x$ such that the Hessian $\mathrm{Hess}_{M_x}^{g,A} f(x) : T_x M_x \to T_x M_x$ can be rewritten as*

$$ (5.16) \qquad \mathrm{Hess}_{M_x}^{g,A} f(x)[\xi] = \tilde{G}_T(x)^{-1} \mathrm{Hess}_{M_x}^{\mathcal{E}} f(x)[\xi] + \tilde{A}_g(\xi, \Pi_x \nabla_{\xi} f), $$

where

(i) $\Pi_x : \mathcal{E} \to T_x M_x$ is the $\mathcal{E}$-orthogonal projection operator on $T_x M_x$;

(ii) $\tilde{G}_T(x)^{-1} : T_x M_x \to T_x M_x$ is the mapping characterized by

$$ \tilde{G}_T(x)^{-1} \xi \in T_x M_x \text{ and } g(\tilde{G}_T(x)^{-1} \xi, \zeta) = \langle \xi, \zeta \rangle_{\mathcal{E}} \text{ for any } \xi, \zeta \in T_x M_x; $$

(iii) $\mathrm{Hess}_{M_x}^{\mathcal{E}}$ is the Riemannian Hessian with respect to the Euclidean metric $g \equiv g^{\mathcal{E}}$ on $M_x$. It reads

$$ \begin{aligned} \mathrm{Hess}_{M_x}^{\mathcal{E}}[\xi] &= \Pi_x \nabla_{\xi}^2 \mathcal{L}(x, \lambda^*(x)) \Pi_x \xi \\ &= \Pi_x [d_{\xi} \nabla_{\xi} f(x) + D_{\xi} \Pi_x ((\mathrm{Id}_{\xi} - \Pi_x) \nabla_{\xi} f(x))], \quad \forall \xi \in T_x M_x, \end{aligned} $$

where $\lambda^*(x) := (Dc(x)Dc(x)^{*\mathcal{E}})^{-1}Dc(x)\nabla_{\mathcal{E}} f(x)$ is the least-squares multiplier.

*Proof.* Due to (5.14), the Hessian $\mathrm{Hess}_{M_x}^{g,A} f(x) : T_x M_x \to T_x M_x$ with respect to the connection $\nabla^A$ and the metric $g$ is the operator defined by

$$ \mathrm{Hess}_{M_x}^{g,A} f(x)[\xi] := \nabla_{\xi}^A \mathrm{grad}_{M_x}^g f = \mathrm{Proj}_x(D_{\xi} \mathrm{grad}_{M_x}^g f) + A(\xi, \mathrm{grad}_{M_x}^g f), \quad \forall \xi \in T_x M_x, $$
where for the moment, we don't specify the tangent projection operator $\text{Proj}_x$. Since $\text{grad}_{\mathcal{M}}^g f = \tilde{G}_T(x)^{-1}\text{Proj}_x^{*,\mathcal{E}}\nabla_{\mathcal{E}}f$, we find

$$ (5.17) \qquad \begin{aligned} \text{Hess}_{\mathcal{M}_x}^{g,A} f(x)[\xi] &= \text{Proj}_x \left( d_{\xi}(\tilde{G}_T(x)^{-1}\text{Proj}_x^{*,\mathcal{E}}\nabla_{\mathcal{E}}f) \right) + A(\xi, \tilde{G}_T(x)^{-1}\text{Proj}_x^{*,\mathcal{E}}\nabla_{\mathcal{E}}f) \\ &= \text{Proj}_x(\tilde{G}_T(x)^{-1}d_{\xi}[\text{Proj}_x^{*,\mathcal{E}}(\nabla_{\mathcal{E}}f)]) \\ &\quad + \text{Proj}_x \left( D_{\xi}(\tilde{G}_T(x)^{-1}\text{Proj}_x^{*,\mathcal{E}})\text{Proj}_x^{*,\mathcal{E}}\nabla_{\mathcal{E}}f \right) + A(\xi, \tilde{G}_T(x)^{-1}\text{Proj}_x^{*,\mathcal{E}}\nabla_{\mathcal{E}}f), \end{aligned} $$

where we note that $\tilde{G}_T(x)^{-1}\text{Proj}_x^{*,\mathcal{E}} = G(x)^{-1}\text{Proj}_x^{*,\mathcal{E}}$ is differentiable as a mapping $\mathcal{E} \to \mathcal{E}$ since $G(x)^{-1}$ and $\text{Proj}_x^{*,\mathcal{E}}$ are differentiable. Let us then recall the Weingarten identities for the dual projector $\text{Proj}_x^{*,\mathcal{E}}$ (Feppon and Lermusiaux, 2019, Proposition 4):

$$ \text{D}_{\eta}\text{Proj}_{x}^{*,\mathcal{E}}(\xi) = (\text{Proj}_{x}^{\perp})^{*,\mathcal{E}}(\text{D}_{\eta}\xi) \text{ for any } \eta \in \text{T}_x\mathcal{M}_x, \xi \in \mathcal{C}^{\infty}(\mathcal{M}_x, \text{N}_x\mathcal{M}_x^{\perp,\mathcal{E}}), $$

$$ \text{D}_{\eta}\text{Proj}_{x}^{*,\mathcal{E}}(\zeta) = -\text{Proj}_{x}^{*,\mathcal{E}}(\text{D}_{\eta}\zeta) \text{ for any } \eta \in \text{T}_x\mathcal{M}_x, \zeta \in \mathcal{C}^{\infty}(\mathcal{M}_x, \text{T}_x\mathcal{M}_x^{\perp,\mathcal{E}}). $$

These identities can be obtained by differentiating the identities $\text{Proj}_x^{*,\mathcal{E}}\xi = \xi$ and $\text{Proj}_x^{*,\mathcal{E}}\zeta = 0$ with respect to $\eta \in \text{T}_x\mathcal{M}_x$ for any smooth vector field $\xi$ with values in $\text{N}_x\mathcal{M}_x^{\perp,\mathcal{E}}$ and any smooth vector field $\zeta$ with values in $\text{T}_x\mathcal{M}_x^{\perp,\mathcal{E}}$. Using Range$(\tilde{G}_T(x)^{-1}) = \text{T}_x\mathcal{M}_x$ and the first Weingarten identity, we obtain

$$ (5.18) \qquad \begin{aligned} \text{Proj}_x(\tilde{G}_T(x)^{-1}d_{\xi}[\text{Proj}_x^{*,\mathcal{E}}\nabla_{\mathcal{E}}f(x)]) \\ &= \tilde{G}_T(x)^{-1}\text{Proj}_x^{*,\mathcal{E}}\text{D}_{\xi}\text{Proj}_x^{*,\mathcal{E}}\nabla_{\mathcal{E}}f(x) + \tilde{G}_T(x)^{-1}\text{Proj}_x^{*,\mathcal{E}}(\text{D}_{\xi}\nabla_{\mathcal{E}}f(x)). \end{aligned} $$

Introducing $\tilde{A}_g(\xi, \zeta) := \text{Proj}_x(\text{D}_{\xi}(\tilde{G}_T(x)^{-1}\text{Proj}_x^{*,\mathcal{E}})\zeta) + A(\xi, \tilde{G}_T(x)^{-1}\zeta)$, (5.17) and (5.18) lead to

$$ \text{Hess}_{\mathcal{M}_x}^{g,A} f(x)[\xi] = \tilde{G}_T(x)^{-1} \text{Proj}_x^{*,\mathcal{E}} \left( \text{D}_{\xi} \nabla_{\mathcal{E}} f(x) + \text{D}_{\xi} \text{Proj}_x^{*,\mathcal{E}} \left( (\text{Proj}_x^{\perp})^{*,\mathcal{E}} \nabla_{\mathcal{E}} f(x) \right) \right) + \tilde{A}_g(\xi, \Pi_x \nabla_{\mathcal{E}} f). $$

The projection perator $\text{Proj}_x$ can be chosen freely in the definition of the covariant derivative by changing the operator $A$ according to lemma 5.3 and when writing $\tilde{G}_T(x)^{-1}\text{Proj}_x^{*,\mathcal{E}}$, according to remark 4.3. The result follows by choosing $\text{Proj}_x = \Pi_x$ to be the $\mathcal{E}$-orthogonal projector $\Pi_x$ on $\text{T}_x\mathcal{M}_x$, which is self-adjoint. The point (iii) can be found in Absil et al. (2009). $\square$

The term in (5.16) featuring $A$ vanishes at critical points $x^*$ satisfying $\Pi_x\nabla_{\mathcal{E}}f(x) = 0$. This implies that the tangent Newton step in (5.13) does not significantly depend on the affine connection $\nabla^A$ or on the tangent metric $G_T(x)$.

**Lemma 5.5.** In the vicinity of a critical point $x^*$ with $\text{grad}_{\mathcal{M}_x^*}^g f(x^*) = 0$ and $\text{Hess}_{\mathcal{M}_x^*}^{\mathcal{E}} f(x^*)$ invertible, all Riemannian Newton steps agree up to a quadratically vanishing error:

$$ (5.19) \qquad -\text{Hess}_{\mathcal{M}_x}^{g,A} f(x)^{-1}\text{grad}_{\mathcal{M}_x}^g f(x) = -(\text{Hess}_{\mathcal{M}_x}^{\mathcal{E}} f(x))^{-1}\Pi_x\nabla_{\mathcal{E}}f(x) + O(||\Pi_x\nabla_{\mathcal{E}}f(x)||_{\mathcal{E}}^2). $$

*Proof.* Let $d_T(x) = -(\text{Hess}_{\mathcal{M}_x}^{g,A} f(x))^{-1}\text{grad}_{\mathcal{M}_x}^g f(x)$. It holds $\text{Hess}_{\mathcal{M}_x}^{g,A} f(x)d_T(x) = -\text{grad}_{\mathcal{M}_x}^g f(x) = -\tilde{G}_T(x)\text{Proj}_x^{*,\mathcal{E}}\nabla_{\mathcal{E}}f(x)$. Using (5.16), left multiplying by $\Pi_x\tilde{G}_T(x)$ and using $\Pi_x\text{Proj}_x^* = \Pi_x$, we find

$$ (\text{Hess}_{\mathcal{M}_x}^{\mathcal{E}} + O(\Pi_x\nabla_{\mathcal{E}}f(x)))d_T(x) = -\Pi_x\nabla_{\mathcal{E}}f(x). $$

The result follows by a standard perturbation analysis. $\square$

We are now in position to prove that Riemannian-Newton landing iterations (5.13) converge quadratically when choosing the normal space $\text{N}_x\mathcal{M}_x = (\nabla_{\mathcal{E}}^2\mathcal{L}(x, \lambda^*(x)))^{\perp,\mathcal{E}}$, characterizing the normal step $d_N(x) = -\text{Dc}(x)^{\dagger,g}c(x) = -\text{Proj}_x^{\perp}\text{Dc}(x)^{\dagger,\mathcal{E}}c(x)$. To make the proof more readable, we don't explicit here the various Lipschitz constants involved in the convergence estimates, absorbing them in the $O(\cdot)$ notation.

**Proposition 5.6.** Assume that $c$ and $f$ are continuous twice differentiable mappings on $\mathcal{E}$. Let $x^* \in \mathcal{M}$ be a local minimizer of (P) with Riemannian Hessian $\text{Hess}_{\mathcal{M}_x^*}^{\mathcal{E}} f(x^*)$ positive definite at $x^*$. If the metric $g$ is set such that the normal space is $\text{N}_{x_k}\mathcal{M}_{x_k} = (\nabla_{\mathcal{E}}^2\mathcal{L}(x_k, \lambda^*(x_k)))^{\perp,\mathcal{E}}$, then, for any $x_0 \in \mathcal{E}$ sufficiently close to $x^*$, the Riemannian-Newton landing iterates (5.13) are well-defined and converge quadratically to $x^*$: there exists a constant $\gamma > 0$ independent of $k \in \mathbb{N}$ such that

$$ ||x_{k+1} - x^*||_{\mathcal{E}} \le \gamma ||x_k - x^*||_{\mathcal{E}}^2, \quad \forall k \in \mathbb{N}. $$
Proof. Due to (5.13) and (5.19).

$$
\begin{align}
(5.20) \qquad & x_{k+1} - x^* = x_k - x^* - \operatorname{Hess}_{M_{x_k}}^{g,A} f(x_k)^{-1} \operatorname{grad}_{M_{x_k}}^g f(x_k) - \operatorname{Dc}(x_k)^{\dagger,g} c(x_k) \nonumber \\
& = x_k - x^* - \operatorname{Hess}_{M_{x_k}}^{\mathcal{E}} f(x_k)^{-1} \Pi_{x_k} \nabla_{\mathcal{E}} f(x_k) - \operatorname{Dc}(x_k)^{\dagger,g} c(x_k) + \mathcal{O}(\|\Pi_{x_k} \nabla_{\mathcal{E}} f(x_k)\|_{\mathcal{E}}^2). \nonumber
\end{align}
$$

Since $\Pi_{x^*} \nabla_{\mathcal{E}} f(x^*) = \text{grad}_{\mathcal{M}_{x^*}}^{\mathcal{E}} f(x^*) = 0$ and $x \mapsto \Pi_x \nabla_{\mathcal{E}} f(x)$ is a Lipschitz map, it is clear that

$$
(5.21) \qquad O(\|\Pi_{x_k} \nabla_{\mathcal{E}} f(x_k)\|_{\mathcal{E}}^2) = O(\|x_k - x^*\|_{\mathcal{E}}^2).
$$

Then, by considering a Taylor expansion, we have on the one hand,

$$
(5.22) \qquad 0 = c(x^*) = c(x_k) + \mathrm{Dc}(x_k)[x^* - x_k] + \mathcal{O}(\|x_k - x^*\|^2_{\mathcal{E}}),
$$

which implies

$$
\operatorname{Proj}_{x_k}^{\perp} (x^* - x_k) = \operatorname{Dc}(x_k)^{\dagger,g} \operatorname{Dc}(x_k)[x^* - x_k] = -\operatorname{Dc}(x_k)^{\dagger,g} c(x_k) + \mathcal{O}(\|x_k - x^*\|^2_{\mathcal{E}}).
$$

Inserting into (5.20), we obtain

$$
(5.23) \qquad x_{k+1} - x^* = \operatorname{Proj}_{x_k}(x_k - x^*) - \operatorname{Hess}_{M_{x_k}}^{\mathcal{E}} f(x_k)^{-1} \Pi_{x_k} \nabla_{\mathcal{E}} f(x_k) + \mathcal{O}(\|x_k - x^*\|^2_{\mathcal{E}}).
$$

On the other hand, performing a Taylor expansion of $x \mapsto \Pi_x \nabla_{\mathcal{E}} f(x)$ around $x_k$ and multiplying by $\Pi_{x_k}$,

$$
(5.24) \qquad 0 = \Pi_{x_k} \Pi_{x^*} \nabla_{\mathcal{E}} f(x^*) = \Pi_{x_k} \Pi_{x_k} \nabla_{\mathcal{E}} f(x_k) + \Pi_{x_k} \mathcal{D}[\Pi_k \nabla_{\mathcal{E}} f](x_k)[x^* - x_k] + \mathcal{O}(\|x_k - x^*\|^2_{\mathcal{E}}).
$$

Recall that for any $\xi \in \mathcal{E}$ (Absil et al., 2009),

$$
\Pi_x \mathcal{D}_\xi (\Pi \cdot \nabla_{\mathcal{E}} f)(x) = \Pi_x (\mathcal{D}_\xi \Pi_x [\Pi_x^\perp \nabla_{\mathcal{E}} f(x)] + \mathcal{D}_\xi \nabla_{\mathcal{E}} f(x)) = \Pi_x \nabla_{\mathcal{E}}^2 \mathcal{L}(x, \lambda^*(x))[\xi].
$$

Consequently,

$$
\begin{align*}
\Pi_{x_k} \mathcal{D}(\Pi \cdot \nabla_{\mathcal{E}} f)(x_k) [x^* - x_k] &= \Pi_{x_k} \nabla_{\mathcal{E}}^2 \mathcal{L}(x_k, \lambda^*(x_k))[x^* - x_k] \\
&= \Pi_{x_k} \nabla_{\mathcal{E}}^2 \mathcal{L}(x_k, \lambda^*(x_k)) \mathrm{Proj}_{x_k}^\perp (x^* - x_k) + \Pi_{x_k} \nabla_{\mathcal{E}}^2 \mathcal{L}(x_k, \lambda^*(x_k)) \mathrm{Proj}_x^\perp (x^* - x_k) \\
&= \mathrm{Hess}_{M_{x_k}}^{\mathcal{E}} f(x_k) \mathrm{Proj}_{x_k}^\perp (x^* - x_k) + \Pi_{x_k} \nabla_{\mathcal{E}}^2 \mathcal{L}(x_k, \lambda^*(x_k)) \mathrm{Proj}_x^\perp (x^* - x_k) \\
&\quad + \mathcal{O}(\|x_k - x^*\|^2_{\mathcal{E}}).
\end{align*}
$$

Substituting into (5.23) after using (5.24), we infer

$$
(5.25) \qquad x_{k+1} - x^* = \mathrm{Hess}_{\mathcal{M}_{x_k}}^{\mathcal{E}} f(x_k)^{-1} \Pi_{x_k} \nabla_{\mathcal{E}}^2 \mathcal{L}(x_k, \lambda^*(x_k)) \mathrm{Proj}_{x_k}^{\perp} (x^* - x_k) + \mathcal{O}(\|x_k - x^*\|^2_{\mathcal{E}}).
$$

Therefore, if $N_{x_k} \mathcal{M}_{x_k} = (\nabla_{\mathcal{E}}^2 \mathcal{L}(x_k, \lambda^*(x_k)))^{\perp, \mathcal{E}}$, we have $\nabla_{\mathcal{E}}^2 \mathcal{L}(x_k, \lambda_k^*) N_{x_k} \mathcal{M}_{x_k} = \mathrm{T}_{x_k} \mathcal{M}_{x_k}^{\perp, \mathcal{E}}$ and

$$
\Pi_{x_k} \nabla_{\mathcal{E}}^2 \mathcal{L}(x_k, \lambda^*(x_k)) \mathrm{Proj}_{x_k}^{\perp} = 0.
$$

It follows that $x_{k+1} - x^* = O(\|x_k - x^*\|^2_\varepsilon)$, which completes the proof. $\square$

If $N_{x_k} \mathcal{M}_{x_k} \neq (\nabla_{\mathcal{E}}^2 \mathcal{L}(x_k, \lambda^*(x_k)))^{\perp, \mathcal{E}}$, (5.25) shows that we only have a priori $x_{k+1} - x^* = O(\|x_k - x^*\|)$, which prevents quadratic convergence.

6. GLOBAL CONVERGENCE USING ADAPTIVE STEP SIZES

In this section, we design an adaptive step size globalization procedure for the landing method.
Namely, we propose a backtracking linesearch based on an Armijo condition for the decrease of a merit
function. Our approach relies heavily on classical techniques for the globalization of SQP (Nocedal
and Wright, 2006; Curtis et al., 2024), leveraging the fact—highlighted in section 5.2—that SQP is a
particular instance of the Riemannian landing.

We consider the landing update with the normal term written in pseudoinverse form, given by

$$
(6.1) \qquad d(x) = d_T(x) + d_N(x) = -\operatorname{grad}_{M_x}^g f(x) - \operatorname{Dc}(x)^{\dagger,g} H(x) c(x)
$$

for some metric $g$ on $\mathcal{E}$ and $H: \mathcal{F} \to \mathcal{F}$, as described in section 4.2.
## 6.1. Assumptions

We introduce a set of generic assumptions on the manifold and Lipschitz constants that will arise in our analysis. The assumptions are standard in the SQP literature (e.g., Berahas et al. (2021); Curtis et al. (2024)). We first assume that there exists a constraint threshold $R > 0$ such that the following compactness assumption holds.

**A1.** There exists a threshold $R$ such that the set $\mathcal{R} = \{x \in \mathcal{E} : \|c(x)\|_{\mathcal{F}} \le R\}$ is compact, contains all the iterates $(x_k)_{k \in \mathbb{N}}$ and trial points, and $f$ is bounded from below on $\mathcal{R}$:

$$ \inf_{x \in \mathcal{R}} f(x) > -\infty. $$

We then require $f$ and $c$ to be Lipschitz continuous on $\mathcal{R}$.

**A2.** $f$ and $c$ are continuously differentiable on $\mathcal{D}$ with Lipschitz derivatives. In particular, there exist constants $||\nabla_{\mathcal{E}}f||_{\infty}, ||Dc||_{\infty}, L_f, L_c > 0$ such that

$$ ||\nabla_{\mathcal{E}}f(x)||_{\mathcal{E}} \le ||\nabla_{\mathcal{E}}f||_{\infty}, \quad \sup_{v \in \mathcal{E} \setminus \{0\}} \frac{||Dc(x)[v]||_{\mathcal{F}}}{||v||_{\mathcal{E}}} \le ||Dc||_{\infty} \quad \forall x \in \mathcal{R}. $$

$$ ||\nabla_{\mathcal{E}}f(x) - \nabla_{\mathcal{E}}f(y)||_{\mathcal{E}} \le L_f ||x - y||_{\mathcal{E}} \quad \forall x, y \in \mathcal{R}, $$

$$ ||Dc(x) - Dc(y)||_{\mathcal{E} \to \mathcal{F}} \le L_c ||x - y||_{\mathcal{E}} \quad \forall x, y \in \mathcal{R}, $$

where $||A||_{\mathcal{E} \to \mathcal{F}} := \sup_{v \in \mathcal{E} \setminus \{0\}} ||Av||_{\mathcal{F}}/||v||$ denotes the operator norm for linear operators $A : \mathcal{E} \to \mathcal{F}$.

Third, we assume the classical linear independence constraint qualification (LICQ) in a neighborhood of the feasible set $\mathcal{M}$:

**A3.** There exists a positive constant $\underline{\sigma}$ such that

$$ \inf_{x \in \mathcal{R}} \sigma_{\min}(Dc(x)) = \inf_{x \in \mathcal{R}} \sigma_m(Dc(x)) \ge \underline{\sigma} > 0, $$

where $\sigma_k(A)$ and $\sigma_{\min}(A)$ denote the $k$th and the smallest (Euclidean) singular value of a linear map $A$, respectively. Equivalently, the Euclidean pseudoinverse of $Dc(x)$ is bounded:

$$ \sup_{w \in \mathcal{F} \setminus \{0\}} \frac{||Dc(x)^{\dagger, \mathcal{E}} w||_{\mathcal{E}}}{||w||_{\mathcal{F}}} \le ||Dc^{\dagger, \mathcal{E}}||_{\infty}, \quad \forall x \in \mathcal{R}. $$

These bounds also imply the uniform boundedness of Riemannian gradients $\nabla_g f(x)$, adjoints $Dc^{*,g}(x)$ and pseudoinverses $Dc^{\dagger,g}$ by assuming uniform ellipticity of the metric.

**A4.** There exist uniform constants $\underline{g} > 0$ and $\bar{g} > 0$ such that

$$ \forall x \in \mathcal{R}, \forall \xi \in \mathcal{E}, \underline{g} \|\xi\|_{\mathcal{E}}^2 \le g(\xi, \xi) = \langle G(x)\xi, \xi \rangle_{\mathcal{E}} \le \bar{g} \|\xi\|_{\mathcal{E}}^2. $$

This assumption implies the following bound which will be useful in the proof of proposition 6.4.

**Lemma 6.1.** Assuming A3, the g- pseudo inverse can be bounded in terms of the Euclidean pseudoinverse:

$$ \sup_{w \in \mathcal{F} \setminus \{0\}} \frac{||Dc(x)^{\dagger,g} w||_{\mathcal{E}}}{||w||_{\mathcal{F}}} \le \sqrt{\frac{\bar{g}}{\underline{g}}} \sup_{w \in \mathcal{F} \setminus \{0\}} \frac{||Dc(x)^{\dagger, \mathcal{E}} w||_{\mathcal{E}}}{||w||_{\mathcal{F}}} \le \sqrt{\frac{\bar{g}}{\underline{g}}} ||Dc^{\dagger, \mathcal{E}}||_{\infty}, \quad \forall x \in \mathcal{R}. $$

*Proof.* This follows from the fact that $Dc(x)^{\dagger,g}w$ is the vector $\xi \in \mathcal{E}$ of smallest g-norm satisfying $Dc(x)\xi = w$. Since $Dc(x)^{\dagger, \mathcal{E}}w$ is another vector satisfying this property, this implies that

$$ \sqrt{\underline{g}} ||Dc(x)^{\dagger,g} w||_{\mathcal{E}} \le ||Dc(x)^{\dagger,g} w||_{g} \le ||Dc(x)^{\dagger, \mathcal{E}} w||_{g} \le \sqrt{\bar{g}} ||Dc(x)^{\dagger, \mathcal{E}} w||_{\mathcal{E}}. $$

□

Finally, we also assume the uniform ellipticity of the operator field $H(x)$ involved in the normal space direction (4.13).

**A5.** There exist uniform constants $\lambda_{\max}(H) > 0$ and $\lambda_{\min}(H) > 0$ such that

$$ \lambda_{\min}(H) ||\zeta||_{\mathcal{F}}^2 \le \langle \zeta, H(x)\zeta \rangle_{\mathcal{F}} \le \lambda_{\max}(H) ||\zeta||_{\mathcal{F}}^2, \quad \text{for all } x \in \mathcal{R} \text{ and } \zeta \in \mathcal{F}. $$

**Remark 6.1.** Assumption A5 holds in both cases $H(x) = \text{Id}_{\mathcal{F}}$ and $H(x) = Dc(x)Dc(x)^{*,g}$.
## 6.2. Line search procedure

The line search procedure is based on decreasing the $l_2$-merit function

$$ (6.2) \qquad \phi_{\mu}(x) = f(x) + \mu \|c(x)\|_{\mathcal{F}} $$

for a large enough penalty parameter $\mu$. We perform a backtracking line-search on the $l_2$-merit $\phi_{\mu}(x)$ with Armijo parameters $\eta \in (0, 1)$ and $\beta \in (0, 1)$ as described in algorithm 1.

**Algorithm 1 Riemannian Landing Method**

**Require:** choose $\eta \in (0, \frac{1}{2})$, $\beta \in (0, 1)$; $\rho \in (0, \lambda_{\min}(H)/2)$; initial $x_0 \in \mathcal{E}$

1: $\mu_0 = 1$

2: $k = 1$

3: **repeat**

4: Compute $d_k = d(x_k)$, the landing direction (6.1) at $x_k$

5: Set $\mu_k = \mu_{k-1}$ if $c(x_k) = 0$ and

$$ (6.3) \qquad \mu_k = \max \left\{ \mu_{k-1}, \frac{Df(x_k)d_N(x_k)}{\rho \|c(x_k)\|_{\mathcal{F}}} \right\} \quad \text{if } c(x_k) \neq 0. $$

6: Set $\alpha_k \leftarrow 1$

7: **while** $\phi_{\mu_k}(x_k + \alpha_k d_k) > \phi_{\mu_k}(x_k) + \eta \alpha_k D\phi_{\mu_k}(x_k) d_k$ **do**

8: $\alpha_k \leftarrow \beta \alpha_k$

▷ backtracking

9: **end while**

10: $x_{k+1} \leftarrow x_k + \alpha_k d_k$

11: **until** convergence

The following lemma gives the directional derivative of the merit function $\phi_{\mu}$.

**Lemma 6.2.** If $x \in \mathcal{E}$ is such that $\|c(x)\| \neq 0$, the merit function $\phi_{\mu}$ is differentiable at $x$ and

$$ (6.4) \qquad D\phi_{\mu}(x)d = Df(x)d + \mu \left\langle \frac{c(x)}{\|c(x)\|}, Dc(x)d \right\rangle_{\mathcal{F}}. $$

For $c(x) = 0$, the merit function is not differentiable at $x$; however, it admits the following directional derivative in any direction $d \in \mathcal{E}$:

$$ (6.5) \qquad D\phi_{\mu}(x)d := \lim_{t \to 0} \frac{\phi_{\mu}(x + td) - \phi_{\mu}(x)}{t} = Df(x)d + \mu \|Dc(x)d\|_{\mathcal{F}}. $$

In what follows, for $x \in \mathcal{R}$, we denote by

$$ (6.6) \qquad d(x) = d_T(x) + d_N(x) $$

the landing step with $d_T(x) = -\mathrm{grad}_{\mathcal{M}_x}^g f(x)$ and $d_N(x) = -\mathrm{Dc}(x)^{\dagger,g} H(x) c(x)$.

The following result shows that if the penalty parameter is large enough, the landing direction is a descent direction for the merit function $\phi_{\mu}$.

**Proposition 6.3** (Sufficient decrease). Let $x \in \mathcal{R}$. If $c(x) = 0$, or if $c(x) \neq 0$ and

$$ (6.7) \qquad \mu \ge \frac{Df(x)d_N(x)}{\rho \|c(x)\|_{\mathcal{F}}}, $$

with $\rho \le \lambda_{\min}(H(x))/2$, then the directional derivative of $\phi_{\mu}(x)$ in the direction $d(x)$ satisfies

$$ (6.8) \qquad D\phi_{\mu}(x)[d(x)] \le -\|\mathrm{grad}_{\mathcal{M}_x}^g f(x)\|_g^2 - \rho\mu\|c(x)\|_{\mathcal{F}}. $$

*Proof.* For $\rho \le \lambda_{\min}(H(x))/2$, we have

$$ (6.9) \qquad \frac{Df(x)d_N(x) \|c(x)\|_{\mathcal{F}}}{\langle c(x), H(x)c(x) \rangle_{\mathcal{F}} - \rho \|c(x)\|_{\mathcal{F}}^2} \le \frac{Df(x)d_N(x) \|c(x)\|_{\mathcal{F}}}{(\lambda_{\min}(H(x)) - \rho) \|c(x)\|_{\mathcal{F}}^2} \le \frac{Df(x)d_N(x)}{\rho \|c(x)\|_{\mathcal{F}}} \le \mu. $$
For $||c(x)|| \neq 0$, (6.4) and $Dc(x)d(x) = -H(x)c(x)$ imply

$$
\begin{align*}
D\phi_{\mu}(x)d(x) &= Df(x)d_{T}(x) + Df(x)d_{N}(x) + \mu \left\langle \frac{c(x)}{\|c(x)\|_{\mathcal{F}}}, Dc(x)d(x) \right\rangle_{\mathcal{F}} \\
&= g(\nabla_{g}f(x), -\mathrm{grad}_{M_{x}}^{g}f(x)) + Df(x)d_{N}(x) - \mu \frac{\langle c(x), H(x)c(x) \rangle_{\mathcal{F}}}{\|c(x)\|_{\mathcal{F}}} \\
&\leq -\|\mathrm{grad}^{g}f(x)\|_{g}^{2} + \frac{\mu}{\|c(x)\|_{\mathcal{F}}} (\langle c(x), H(x)c(x) \rangle_{\mathcal{F}} - \rho \|c(x)\|_{\mathcal{F}}^{2}) - \mu \frac{\langle c(x), H(x)c(x) \rangle_{\mathcal{F}}}{\|c(x)\|_{\mathcal{F}}} \\
&= -\|\mathrm{grad}^{g}f(x)\|_{g}^{2} - \mu\rho \|c(x)\|_{\mathcal{F}},
\end{align*}
$$

where we have used (6.9) in the third line. For $x \in \mathcal{M}$, where $c(x) = 0$, the landing direction $d$ only has a tangent component ($d_N(x) = 0$) and the directional derivative (6.5) becomes

$$D\phi_{\mu}(x)d(x) = Df(x)d_{T}(x) = -\|\mathrm{grad}^{g}f(x)\|_{g}^{2}.$$

Thus (6.8) holds for all $x \in \mathcal{R}$.

The following result shows that the penalty parameter $\mu$ remains bounded in the region of interest $\mathcal{R}$.

**Proposition 6.4** (Upper bound on the penalty parameter). Under the above assumptions, the increasing merit parameter sequence $(\mu_k)_{k \ge 0}$ defined by (6.3) in algorithm 1 remains bounded:

$$(6.10) \qquad \mu(x) := \frac{Df(x)d_N(x)}{\rho \|c(x)\|_{\mathcal{F}}} \le \bar{\mu} := \frac{1}{\rho} \sqrt{\frac{\bar{g}}{g}} \|\nabla_{\mathcal{E}}f\|_{\infty} \|Dc^{\dagger, \mathcal{E}}\|_{\infty} \lambda_{\max}(H), \quad \forall x \in \mathcal{R}.$$

In particular, the sequence $(\mu_k)_{k \in \mathbb{N}}$ is convergent.

*Proof.* Recall that $\mu(x)$. Using the Cauchy-Schwartz inequality and lemma 6.1, we find

$$
\begin{align*}
|Df(x)d_N(x)| &= \langle \nabla_{\mathcal{E}}f(x), -Dc(x)^{\dagger,g}H(x)c(x) \rangle_{\mathcal{E}} \\
&\le \sqrt{\frac{\bar{g}}{g}} \|\nabla_{\mathcal{E}}f(x)\|_{\mathcal{E}} \|Dc^{\dagger,\mathcal{E}}\|_{\infty} \|H(x)c(x)\|_{\mathcal{E}} \\
&\le \sqrt{\frac{\bar{g}}{g}} \|\nabla_{\mathcal{E}}f\|_{\infty} \|Dc^{\dagger,\mathcal{E}}\|_{\infty} \lambda_{\max}(H) \|c(x)\|_{\mathcal{E}}.
\end{align*}
$$

In what follows, we use the following inequality, actually valid for any Euclidean norm.

**Lemma 6.5.** Let $v, w \in \mathcal{F}$ with $v \neq 0$. The following inequality holds:

$$(6.11) \qquad \|v+w\|_{\mathcal{F}} \le \|v\|_{\mathcal{F}} + \left\langle \frac{v}{\|v\|_{\mathcal{F}}}, w \right\rangle_{\mathcal{F}} + \frac{\|w\|_{\mathcal{F}}^2}{2\|v\|_{\mathcal{F}}}.$$

*Proof.* We have

$$(6.12) \qquad \|v+w\|_{\mathcal{F}}^2 = \|v\|_{\mathcal{F}}^2 + 2 \langle v, w \rangle_{\mathcal{F}} + \|w\|_{\mathcal{F}}^2 = \|v\|_{\mathcal{F}}^2 (1+t).$$

with $t := (2 \langle v, w \rangle + \|w\|_{\mathcal{F}}^2) / \|v\|_{\mathcal{F}}^2$. Since $1+t = \|v+w\|_{\mathcal{F}}^2 / \|v\|_{\mathcal{F}}^2 \ge 0$, it is clear that $t \ge -1$. Recalling the elementary inequality $\sqrt{1+t} \le 1+t/2$, valid for $t \ge -1$, this implies

$$(6.13) \qquad \|v+w\|_{\mathcal{F}} = \|v\|_{\mathcal{F}} \sqrt{1+t} \le \|v\|_{\mathcal{F}} \left(1+\frac{t}{2}\right) = \|v\|_{\mathcal{F}} \left(1+\frac{\langle v,w \rangle}{\|v\|_{\mathcal{F}}^2}+\frac{\|w\|_{\mathcal{F}}^2}{2\|v\|_{\mathcal{F}}^2}\right).$$

The following lemma gives an upper bound on the merit function along the landing direction, which we use below to show finite termination of the linesearch.

**Lemma 6.6** (Quadratic upper bound for the merit function). For any $x \in \mathcal{R}$ and $\alpha \ge 0$ such that $x + \alpha d(x) \in \mathcal{R}$ with $d(x)$ the landing step (6.6), it holds that

$$(6.14) \qquad \phi_{\mu}(x + \alpha d(x)) \le \phi_{\mu}(x) + \alpha D\phi_{\mu}(x)d(x) + \frac{\alpha^2}{2}(L_f + \mu L_c) \|d(x)\|_{\mathcal{E}}^2 + \frac{\alpha^2}{2}\mu \lambda_{\max}(H)^2 \|c(x)\|_{\mathcal{F}},$$

where $D\phi_{\mu}(x)d$ denotes the one-sided directional derivative of $\phi_{\mu}$ at $x$ along $d$, given by Lemma 6.2.
Proof. By the Lipschitz continuity of $\nabla_{\mathcal{E}}f$ (assumption A2), we can write

$$ (6.15) \qquad f(x + \alpha d(x)) \le f(x) + \alpha \mathrm{D}f(x)d(x) + \frac{1}{2} L_f \alpha^2 \|d(x)\|_{\mathcal{E}}^2. $$

By the Lipschitz continuity of $Dc$ and $Dc(x)d(x) = Dc(x)d_N(x) = -H(x)c(x)$, we also have

$$ (6.16) \qquad \begin{aligned} c(x + \alpha d(x)) &= c(x) + \alpha Dc(x)d(x) + r(\alpha), \\ &= (\mathrm{Id}_{\mathcal{F}} - \alpha H(x))c(x) + r(\alpha) \end{aligned} $$

with $\|r(\alpha)\|_{\mathcal{F}} \le \frac{1}{2} L_c \alpha^2 \|d(x)\|_{\mathcal{E}}^2$.

If $c(x) = 0$, it follows that $\|c(x + \alpha d)\|_{\mathcal{F}} \le \|r(\alpha)\|_{\mathcal{F}}$, which imply

$$ (6.17) \qquad \begin{aligned} \phi_{\mu}(x + \alpha d(x)) &= f(x + \alpha d(x)) + \mu \|c(x + \alpha d(x))\|_{\mathcal{F}} \\ &\le f(x) + \alpha \mathrm{D}f(x)d(x) + \frac{1}{2} L_f \alpha^2 \|d(x)\|_{\mathcal{E}}^2 + \mu \frac{1}{2} L_c \alpha^2 \|d(x)\|_{\mathcal{E}}^2. \end{aligned} $$

Moreover, (6.5) together with $Dc(x)d(x) = 0$ entail $\mathrm{D}f(x)d(x) = \mathrm{D}\phi_{\mu}(x)d(x)$. Therefore (6.14) holds if $c(x) = 0$.

We now consider the case $c(x) \ne 0$. Majoring the first line of section 6.2 using (6.11) with $v = c(x)$ and $w = \alpha Dc(x)d(x)$ gives

$$ (6.18) \qquad \begin{aligned} \|c(x+\alpha d(x))\|_{\mathcal{F}} &\le \|c(x)\|_{\mathcal{F}} + \alpha \left\langle \frac{c(x)}{\|c(x)\|_{\mathcal{F}}}, \mathrm{D}c(x)d(x) \right\rangle_{\mathcal{F}} + \frac{\alpha^2}{2} \frac{\|\mathrm{D}c(x)d(x)\|_{\mathcal{F}}^2}{\|c(x)\|_{\mathcal{F}}} + \|r(\alpha)\|_{\mathcal{F}} \\ &\le \|c(x)\|_{\mathcal{F}} + \alpha \left\langle \frac{c(x)}{\|c(x)\|_{\mathcal{F}}}, \mathrm{D}c(x)d(x) \right\rangle_{\mathcal{F}} + \frac{\alpha^2}{2} \lambda_{\max}(H)^2 \|c(x)\|_{\mathcal{F}} + \frac{1}{2} L_c \alpha^2 \|d\|_{\mathcal{E}}^2. \end{aligned} $$

Combining (6.18), (6.15) and (6.4) yields (6.14). $\square$

The following result shows that the linesearch terminates in a finite number of steps, and gives an explicit lower bound on the stepsize.

**Proposition 6.7** (Lower bound on the step size). The sequence of step sizes $(\alpha_k)_{k \in \mathbb{N}}$ generated by algorithm 1 are bounded from below:

$$ (6.19) \qquad \alpha_k \ge \underline{\alpha} := \frac{2\beta(1-\eta)}{(L_f + \bar{\mu}L_c)} \min\left\{1, \frac{\rho}{c_1}\right\} > 0 \quad \forall k \in \mathbb{N}, $$

where $c_1$ is the constant

$$ (6.20) \qquad c_1 := \left( \bar{g} \|\mathrm{D}c^{\dagger, \mathcal{E}}\|_{\infty}^2 R + \frac{g\bar{\mu}}{(L_f + L_c)} \right) \lambda_{\max}(H)^2. $$

Proof. Owing to (6.14), the Armijo condition

$$ (6.21) \qquad \phi_{\mu}(x + \alpha d(x)) \le \phi_{\mu}(x) + \eta \alpha D\phi_{\mu}(x)d(x) $$

is satisfied whenever

$$ (6.22) \qquad \alpha D\phi_{\mu}(x)d(x) + \frac{\alpha^2}{2}(L_f + \mu L_c)\|d\|_{\mathcal{E}}^2 + \frac{\alpha^2}{2}\mu\lambda_{\max}(H)^2 \|c(x)\|_{\mathcal{F}} \le \eta\alpha D\phi_{\mu}(x)d(x), $$

or, equivalently,

$$ (6.23) \qquad \alpha \le \alpha_{\min} \text{ with } \alpha_{\min} := \frac{2(1-\eta)|\mathrm{D}\phi_{\mu}(x)d(x)|}{(L_f + \mu L_c) \|d(x)\|_{\mathcal{E}}^2 + \mu \lambda_{\max}(H)^2 \|c(x)\|_{\mathcal{F}}}. $$

This means that, at iteration $k$, the Armijo line search terminates with $\alpha_k > 0$ satisfying either $\alpha_k = 1$ or $\alpha_k/\beta \ge \alpha_{\min}$. The sufficient decrease condition (6.8) implies then

$$ (6.24) \qquad \alpha_k \ge \frac{2\beta(1-\eta)}{L_f + \mu_k L_c} \times \frac{\|\mathrm{grad}_{M_{x_k}}^g f(x_k)\|_g^2 + \rho\mu_k \|c(x_k)\|_{\mathcal{F}}}{\|d_k\|_{\mathcal{E}}^2 + \mu_k \lambda_{\max}(H)^2 \|c(x_k)\|_{\mathcal{F}} / (L_f + \mu_k L_c)}. $$

By using A4 and the $g$-orthogonality of the tangent and normal step, we descent direction $d_k$ can be estimated as

$$ \begin{aligned} \|d_k\|_{\mathcal{E}}^2 &\le \frac{1}{g} g(d_k, d_k) = \frac{1}{g} (g(d_T(x_k), d_T(x_k)) + g(d_N(x_k), d_N(x_k))) \\ &\le \frac{1}{g} \|\mathrm{grad}_{M_{x_k}}^g f(x_k)\|_g^2 + \frac{1}{g} \|\mathrm{D}c(x_k)^{\dagger, \mathcal{E}} H(x_k) c(x_k)\|_g^2 \\ &\le \frac{1}{g} \left( \|\mathrm{grad}_{M_{x_k}}^g f(x_k)\|_g^2 + \bar{g} \|\mathrm{D}c^{\dagger, \mathcal{E}}\|_{\infty}^2 \lambda_{\max}(H)^2 \|c(x_k)\|_{\mathcal{F}}^2 \right), \end{aligned} $$
where we use the property that $d_N(x_k)$ is the minimum $g$-norm vector $\xi \in \mathcal{E}$ satisfying $Dc(x_k)\xi = -H(x_k)c(x_k)$ and that $-Dc^{\dagger,\mathcal{E}}(x_k)H(x_k)c(x_k)$ is another such vector. Therefore, for $\|c(x)\|_{\mathcal{F}} \le R$, the denominator of (6.24) is upper bounded as follows

$$ (6.25) \qquad \begin{aligned} & \frac{\|d_k\|_{\mathcal{E}}^2 + \mu_k \lambda_{\max}(H) \|c(x_k)\|_{\mathcal{F}}}{L_f + \mu_k L_c} \\ & \le \frac{1}{g} \left\| \mathrm{grad}_{M_{x_k}}^g f(x_k) \right\|_g^2 + \frac{\bar{g}}{g} \left\| Dc^{\dagger,\mathcal{E}} \right\|_{\infty}^2 \lambda_{\max}(H)^2 \|c(x_k)\|_{\mathcal{F}}^2 + \frac{\mu_k \lambda_{\max}(H)^2}{L_f + \mu_k L_c} \|c(x_k)\|_{\mathcal{F}} \\ & \le \frac{1}{g} \left\| \mathrm{grad}_{M_{x_k}}^g f(x_k) \right\|_g^2 + \left( \frac{\bar{g}}{g} \left\| Dc^{\dagger,\mathcal{E}} \right\|_{\infty}^2 R + \frac{\mu_k}{L_f + \mu_k L_c} \right) \lambda_{\max}(H)^2 \|c(x_k)\|_{\mathcal{F}} \\ & \le \frac{1}{g} \left( \left\| \mathrm{grad}_{M_{x_k}}^g f(x_k) \right\|_g^2 + c_1 \|c(x_k)\|_{\mathcal{F}} \right), \end{aligned} $$

where $c_1$ is defined in (6.20) and the last inequality follows from $1 \le \mu_k \le \bar{\mu}$ for all $k$. Thus we have

$$ (6.26) \qquad \alpha_k \ge \frac{2\beta(1-\eta)}{(L_f + \bar{\mu}L_c)} F(\|\mathrm{grad}_{M_{x_k}}^g f(x_k)\|_g^2, \|c(x_k)\|_{\mathcal{F}}) \quad \text{with } F(t_1, t_2) := \frac{t_1 + \rho t_2}{t_1 + c_1 t_2}. $$

Finally, for $t_1, t_2 \ge 0$ with $t_1 + c_1 t_2 > 0$, it holds $F(t_1, t_2) \ge \min(1, \rho/c_1)$. This implies the bound. $\square$

We are now in a position to prove our main convergence result for the Riemannian landing method. The convergence rate is consistent with an existing result for SQP (Curtis et al., 2024, Thm. 1).

**Theorem 6.8 (Global convergence of the landing iterates).** Under assumptions A1 to A5, the iterates $(x_k)_{k \in \mathbb{N}}$ generated by algorithm 1 satisfy

$$ (6.27) \qquad \|c(x_k)\|_{\mathcal{F}} \to 0 \quad \text{and} \quad \|\mathrm{grad}_{M_x}^g f(x_k)\|_{\mathcal{E}} \to 0 \quad \text{as } k \to +\infty. $$

Additionally, the rate of convergences of these sequences can be estimated as follows:

$$ (6.28) \qquad \min_{k=0,...,K-1} \|\mathrm{grad}_{M_x}^g f(x_k)\|_g \le \frac{C}{\sqrt{K}} $$

and

$$ (6.29) \qquad \min_{k=0,...,K-1} \|c(x_k)\|_{\mathcal{F}} \le \frac{C}{K}. $$

for a constant $C$ independent of $K$.

*Proof.* For $K \ge 0$, the Armijo condition reads

$$ \phi_{\mu_k}(x_{k+1}) \le \phi_{\mu_k}(x_k) + \eta \alpha_k D \phi_{\mu_k}(x_k) d_k, $$

which, thanks to (6.8), implies

$$ f(x_{k+1}) + \mu_k \|c(x_{k+1})\|_{\mathcal{F}} \le f(x_k) + \mu_k \|c(x_k)\|_{\mathcal{F}} - \eta \alpha_k \|\mathrm{grad}_{M_{x_k}}^g f(x_k)\|_g^2 - \eta \alpha_k \rho \mu_k \|c(x_k)\|_{\mathcal{F}}. $$

This can be rewritten as

$$ (6.30) \qquad \begin{aligned} & f(x_k) - f(x_{k+1}) + (\mu_k \|c(x_k)\|_{\mathcal{F}} - \mu_{k+1} \|c(x_{k+1})\|_{\mathcal{F}}) \\ & \ge (\mu_k - \mu_{k+1}) \|c(x_{k+1})\|_{\mathcal{F}} + \eta \alpha_k \|\mathrm{grad}_{M_{x_k}}^g f(x_k)\|_g^2 + \eta \alpha_k \rho \mu_k \|c(x_k)\|_{\mathcal{F}}. \end{aligned} $$

Summing for $0 \le k \le K-1$ and using the lower bound (6.19), we obtain

$$ \begin{aligned} & f(x_0) - f(x_T) + \mu_0 \|c(x_0)\|_{\mathcal{F}} \\ & \ge \mu_T \|c(x_T)\|_{\mathcal{F}} + \sum_{k=1}^{K} (\mu_{k-1} - \mu_k) \|c(x_k)\|_{\mathcal{F}} + \sum_{k=0}^{K-1} \eta \alpha_k \rho \mu_k \|c(x_k)\|_{\mathcal{F}} + \sum_{k=0}^{K-1} \eta \alpha_k \|\mathrm{grad}_{M_{x_k}}^g f(x_k)\|_g^2 \\ & \ge \sum_{k=1}^{K-1} (\eta \alpha \rho + \mu_{k-1} - \mu_k) \|c(x_k)\|_{\mathcal{F}} + \mu_{K-1} \|c(x_{K-1})\|_{\mathcal{F}} + \sum_{k=0}^{K-1} \eta \alpha_k \|\mathrm{grad}_{M_{x_k}}^g f(x_k)\|_g^2. \end{aligned} $$

Since $\mu_{k-1} - \mu_k \to 0$ as $k \to +\infty$, we conclude that the series $\sum_k \|\mathrm{grad}_{M_{x_k}}^g f(x_k)\|_g^2$ and $\sum_k \|c(x_k)\|_{\mathcal{F}}$ are convergent, which implies (6.27) to (6.29). $\square$

Due to the boundedness of the sequence $(x_k)_{k \ge 0}$, there must be an accumulation point, which is necessary a feasible KKT point according to (6.27).
## 7. GEOMETRIC DESIGN OF TANGENT AND NORMAL TERMS FOR MATRIX OPTIMIZATION WITH ORTHOGONALITY CONSTRAINTS

A natural application for the landing algorithm is optimization with orthogonality constraints, such as

$$ (7.1) \qquad \min_{X \in \mathbb{R}^{n \times p}} f(X) \text{ subject to } X^\top X = I_p, $$

where $p \le n$. The feasible set—called the Stiefel manifold—consists of rectangular matrices with orthonormal columns

$$ \mathcal{M} = \mathrm{St}(n, p) := \{X \in \mathbb{R}^{n \times p} | X^\top X = I_p\}. $$

A function that defines the constraints is

$$ (7.2) \qquad c: \mathbb{R}^{n \times p} \to \mathrm{Sym}(p): c(X) = \frac{1}{2}(X^\top X - I_p). $$

Thus, we have an instance of (P) with $\mathcal{E} = \mathbb{R}^{n \times p}$ and $\mathcal{F} = \mathrm{Sym}(p)$, the set of symmetric matrices of size $p$. The spaces $\mathcal{E}$ and $\mathcal{F}$ are equipped with the standard Frobenius inner product $\langle \cdot, \cdot \rangle_{\mathcal{E}}$ and $\langle \cdot, \cdot \rangle_{\mathcal{F}}$. In all what follows, we drop the subscript notation $\cdot_{\mathcal{E}}$ and $\cdot_{\mathcal{F}}$ when writing matrix Frobenius inner products or norms, in particular, for two matrices of same size $X$ and $Y$, $\langle X, Y \rangle = \mathrm{trace}(X^\top Y)$ and $||X||^2 = \langle X, X \rangle$. Let $\mathbb{R}_*^{n \times p}$ denote the matrices size $n \times p$ with full rank. We write $\mathrm{Sym}(p) = \{S \in \mathbb{R}^{p \times p} | S^\top = S\}$ and $\mathrm{Skew}(p) = \{\Omega \in \mathbb{R}^{p \times p} | \Omega^\top = -\Omega\}$, and $\mathrm{sym}(A) = (A + A^\top)/2$ and $\mathrm{skew}(A) = (A - A^\top)/2$.

The squared infeasibility for the problem (7.1) reads

$$ \psi(X) = \frac{1}{2} \|c(X)\|^2 = \frac{1}{4} \|X^\top X - I_p\|^2, $$

where $I_p$ is the identity matrix of $\mathbb{R}^{p \times p}$. The set $\mathcal{D}$ of (3.1) is the set of full rank matrices $X \in \mathbb{R}_*^{n \times p}$ and the layer manifold $\mathcal{M}_X$ is

$$ \mathcal{M}_X = \mathrm{St}_{X^\top X} := \{Y \in \mathbb{R}^{n \times p} : Y^\top Y = X^\top X\}. $$

In this section, we demonstrate how the family of normal spaces $X \to N_X \mathrm{St}_{X^\top X}$ is crucial to lead to efficient computations of the tangent and normal terms of the landing algorithm (1.2). To illustrate this claim, we first consider the family of normal spaces orthogonal to the tangent spaces with respect to the Euclidean Frobenius inner product. Although this choice is natural, computing the corresponding orthogonal projections is expensive because it requires solving Sylvester equations. We then take the reverse approach: we first introduce an alternative family of normal spaces with associated projection operators that can be computed explicitly, then we design tangent and normal metrics that lead to explicit—and thus more exploitable—formulas for the tangent and normal steps for matrix optimization with orthogonality constraints.

The starting point is to recall the characterization of the tangent space to $\mathrm{St}_{X^\top X}$ at $X \in \mathbb{R}_*^{n \times p}$. We emphasize that this set does not depend on the metric.

**Proposition 7.1** (Gao et al. (2022); Goyens et al. (2026)). The tangent space of $\mathrm{St}_{X^\top X}$ at $X$ is the set

$$ (7.3) \qquad \mathrm{T}_X \mathrm{St}_{X^\top X} = \{\xi \in \mathbb{R}^{n \times p} | \xi^\top X + X^\top \xi = 0\} $$

$$ (7.4) \qquad = \{X(X^\top X)^{-1}\Omega + \Delta | \Omega \in \mathrm{Skew}(p), \Delta \in \mathbb{R}^{n \times p} \text{ with } \Delta^\top X = 0\} $$

$$ (7.5) \qquad = \{WX | W \in \mathrm{Skew}(n)\}, $$

with dimension $np - p(p + 1)/2$.

### 7.1. Landing algorithm for the Euclidean metric

We first outline the derivation of the tangent and normal steps (1.3) and (4.13) when $g \equiv g^\mathcal{E}$ is the Euclidean or Frobenius inner product metric:

$$ g^\mathcal{E}(\xi, \zeta) = \langle \xi, \zeta \rangle \quad \forall \xi, \zeta \in \mathbb{R}^{n \times p}. $$

We thus choose the formulation of the landing algorithm (4.12) based on the pseudoinverse formula for the normal step. The tangent and normal vector fields (1.3) and (4.13), with respect to the Euclidean metric, are given by

$$ (7.6) \qquad d_T(X) = -\mathrm{Proj}_{X, \mathcal{E}}[\nabla_{\mathcal{E}} f(X)], $$

$$ (7.7) \qquad d_N(X) = -\mathrm{Dc}(X)^{\dagger, \mathcal{E}}[H(X)[c(X)]] $$
where brackets denote the action of linear operators on matrix spaces.

In what follows, we compute these two directions explicitly. The orthogonal projector $\text{Proj}_{X,\mathcal{E}}$ can be obtained through the characterization (3.5), as is commonly done in the literature on Riemannian optimization Absil et al. (2008). However, since the formula for $d_N(X)$ also requires the computation of the pseudoinverse $\text{Dc}(X)^{\dagger,\mathcal{E}}$, we instead propose to compute it using the explicit relation $\text{Proj}_{X,\mathcal{E}} = \text{Id}_{\mathcal{E}} - \text{Dc}(X)^{\dagger,\mathcal{E}}\text{Dc}(X)$. As it becomes clear below, the main difficulty in computing this projection lies in inverting the symmetric operator $(\text{Dc}(X)\text{Dc}(X)^{*},\mathcal{E})^{-1}$ that appears in the pseudoinverse formula $\text{Dc}(X)^{\dagger,\mathcal{E}} = \text{Dc}(X)^{*},\mathcal{E}(\text{Dc}(X)\text{Dc}(X)^{*},\mathcal{E})^{-1}$. Let us start by evaluating $\text{Dc}(X)$ and $\text{Dc}(X)^{*},\mathcal{E}$.

**Proposition 7.2.** The operator $\text{Dc}(X) : \mathbb{R}^{n \times p} \to \text{Sym}(p)$ and its adjoint $\text{Dc}(X)^{*},\mathcal{E} : \text{Sym}(p) \to \mathbb{R}^{n \times p}$ read:

(i) $\text{Dc}(X)[\xi] = \frac{1}{2}(X^{\top}\xi + \xi^{\top}X) = \text{sym}(X^{\top}\xi)$ for all $\xi \in \mathbb{R}^{n \times p}$,

(ii) $\text{Dc}(X)^{*},\mathcal{E}[S] = XS$ for all $S \in \text{Sym}(p)$.

*Proof.* The point (i) is obvious in view of (7.2). For (ii), we write, for any $S \in \text{Sym}(p)$ and $\xi \in \mathbb{R}^{n \times p}$,

$$ \langle \xi, \text{Dc}(X)^{*},\mathcal{E}[S] \rangle = \langle \text{Dc}(X)[\xi], S \rangle = \langle \text{sym}(X^{\top}\xi), S \rangle = \langle X^{\top}\xi, S \rangle = \langle \xi, XS \rangle. $$

□

We infer the following characterization of the normal space $N_{X}^{\mathcal{E}}\text{St}_{X^{\top}X} = \text{Range}(\text{Dc}(X)^{*},\mathcal{E})$ and of the Euclidean pseudoinverse $\text{Dc}(X)^{\dagger,\mathcal{E}}$.

**Corollary 7.3.** The normal space of $\text{St}_{X^{\top}X}$ at $X \in \mathbb{R}_{*}^{n \times p}$ with respect to the Euclidean metric $g^{\mathcal{E}}$ is
$$ (7.8) \qquad N_{X}^{\mathcal{E}}\text{St}_{X^{\top}X} = \{XS \mid S \in \text{Sym}(p)\}. $$

**Proposition 7.4.** The pseudoinverse $\text{Dc}(X)^{\dagger,\mathcal{E}} : \text{Sym}(p) \to \mathbb{R}^{n \times p}$ is the operator defined by

$$ \text{Dc}(X)^{\dagger,\mathcal{E}}[T] = XS, \quad \forall T \in \text{Sym}(p), $$

where $S$ is the unique solution in $\text{Sym}(p)$ to the Sylvester equation

$$ (7.9) \qquad \frac{1}{2}(X^{\top}XS + SX^{\top}X) = T. $$

*Proof.* From proposition 7.2, it is readily seen that

$$ \text{Dc}(X)\text{Dc}(X)^{*},\mathcal{E}[S] = \frac{1}{2}(X^{\top}XS + SX^{\top}X), $$

from where these results follow easily.

□

Consequently, the computation of the tangent and the normal terms of the landing algorithm require, in full generality, solving Sylvester equations.

**Corollary 7.5.** (i) The tangent term (7.6) of the landing algorithm (1.3) with respect to the Euclidean metric given by

$$ (7.10) \qquad d_T(X) = -\text{Proj}_{X,\mathcal{E}}[\nabla_{\mathcal{E}}(f(X))] = -\nabla_{\mathcal{E}}f(X) + XS, \quad \forall Z \in \mathbb{R}^{n \times p}, $$

where $S$ is the unique solution in $\text{Sym}(p)$ to the Sylvester equation (7.9) with $T = \text{sym}(X^{\top}\nabla_{\mathcal{E}}f(X))$.

(ii) The normal vector field (7.7) is given by $d_N(X) = -XS$, where $S$ is the solution to the Sylvester equation (7.9) with $T = \frac{1}{2}H(X)[X^{\top}X - I_p]$.

*Proof.* These results are obtained by using proposition 7.4 and formula (3.15) for the projection operator $\text{Proj}_{X,\mathcal{E}}$. □

In the general case, we mention that the solution to (7.9) could be calculated from the singular value decomposition of $X$. The following result is classical and may be found e.g. in (Li and Zhou, 2018).

**Lemma 7.6.** Let $T \in \mathbb{R}^{p \times p}$ be an arbitrary matrix. Let $X = \sum_{i=1}^{p} \sigma_i u_i v_i^{\top}$ be the singular value decomposition of a full rank matrix $X \in \mathbb{R}^{n \times p}$ with singular values $\sigma_1 \ge \sigma_2 \ge \dots \ge \sigma_p > 0$, left and right singular vectors $(u_i)_{1 \le i \le p}$ and $(v_i)_{1 \le i \le p}$. The unique solution $S \in \mathbb{R}^{p \times p}$ to the Sylvester equation (7.9) is given by

$$ S = \frac{1}{2} \sum_{1 \le i \le j \le p} \frac{v_i^{\top} T v_j}{\sigma_i^2 + \sigma_j^2} (v_i v_j^{\top} + v_j v_i^{\top}). $$
When $T = Q(X^\top X)$ with $Q$ an analytic function, the solution to the Sylvester equation (7.9) is explicitly given by

$$S = Q(X^\top X)(X^\top X)^{-1}.$$

This observation leads to explicit formulas for the normal field $d_N(X)$ for two particular values of the operator $H(X) : \text{Sym}(p) \to \text{Sym}(p)$.

**Proposition 7.7.** (i) If $H(X) = \text{Id}_F$, then $d_N(X) = -\frac{1}{2}X(I_p - (X^\top X)^{-1})$.

(ii) if $H(X) = \text{Dc}(X)\text{Dc}(X)^*, \varepsilon$, then $d_N(X) = -\nabla_\varepsilon\psi(X) = -X(X^\top X - I_p)$.

*Proof.* For $H(X) = \text{Id}_F$, $d_N(X) = -XS$ with $S$ being the solution to (7.9) with $T = Q(X^\top X)$ with $Q(x) = \frac{1}{2}(x-1)$. The solution is $S = (I_p - (X^\top X)^{-1})/2$, which yields (i). For $H(X) = \text{Dc}(X)\text{Dc}(X)^*, \varepsilon$, (ii) can be found directly by using the identity $d_N(X) = -\nabla_\varepsilon\psi(X)$ with $\psi(X) = \frac{1}{4}\|X^\top X - I_p\|_F^2$. It is instructive, however, to retrieve this result from the result of corollary 7.5. We have in this case that $S$ is the solution to (7.9) the right-hand side $T$ given by

$$ \begin{aligned} T &= \frac{1}{2}\text{Dc}(X)\text{Dc}(X)^*, \varepsilon [X^\top X - I_p] \\ &= \frac{1}{2}[X^\top(XX^\top X - X) + (X^\top XX^\top - X^\top)X] \\ &= (X^\top X)^2 - X^\top X \\ &= Q(X^\top X) \text{ with } Q(x) = x^2 - x. \end{aligned} $$

Hence, the solution to the Sylvester equation (7.9) is $S = (X^\top X - I_p)$, which yields (ii). $\square$

To summarize, the execution of the landing algorithm (1.3) for minimizing a function on the Stiefel manifold using the Euclidean metric $g^\varepsilon$ to define the tangent and the normal terms $d_T(X)$ and $d_N(X)$ requires solving a Sylvester equation for computing the tangent term at every iteration, which is potentially costly for large matrices $X$.

## 7.2. Metric design based on an explicit choice of projection operators

In this section, we adopt the reverse approach: instead of first defining the metric $g$ and then deducing the associated tangent and normal terms $d_T(X)$ and $d_N(X)$ via (4.3) and (4.13), we begin by choosing first an explicit family of projection operators $\text{Proj}_X$ onto the tangent space $T_X\text{St}_{X^\top X}$, or equivalently, a family of normal spaces $X \mapsto N_X\text{St}_{X^\top X}$. This step is outlined in details in section 7.2.1. In the second step, we construct explicit operators

$$ (7.11) \qquad \widetilde{G}_T(X) : T_X\text{St}_{X^\top X} \to N_X\text{St}_{X^\top X}^\perp, \qquad \widetilde{G}_N(X) : N_X\text{St}_{X^\top X} \to T_X\text{St}_{X^\top X}, $$

each admitting explicit inverses. This construction defines the metric a posteriori through (4.6):

$$ (7.12) \qquad g(\xi, \zeta) = \langle \widetilde{G}_T(X)\text{Proj}_X[\xi], \zeta \rangle + \langle \widetilde{G}_N(X)\text{Proj}_X^\perp[\xi], \zeta \rangle. $$

Moreover, these explicit inverses provide closed-form expressions for the tangent and normal terms of the landing algorithm via (4.10) and (4.17).

We examine two possible choices for $\widetilde{G}_T$ and $\widetilde{G}_N$. In section 7.2.2, we consider a 'canonical' choice for these mappings, which leads to the definition of a new metric $g$ on $\mathbb{R}_*^{n \times p}$ and new formulas for the landing algorithm for the optimization problem (7.1). In section 7.2.3, we consider instead the family of $\beta$-metrics $g^\beta$ introduced in Goyens et al. (2026), which includes for $\beta = \frac{1}{2}$ the pulled-back canonical metric for the layered Stiefel manifold studied in Gao et al. (2022). We show that this family of metric is associated to the proposed family of normal spaces $X \mapsto N_X\text{St}_{X^\top X}$ considered in the first step, with a specific choice of $\widetilde{G}_T$ and $\widetilde{G}_N$. These operators turn to admit explicit inverses, thereby highlighting why closed-form formulas for $d_T(X)$ and $d_N(X)$ are also available in this case.

### 7.2.1. Choice of oblique projection operators and associated normal spaces

Recalling that any tangent vector $\xi \in T_X\text{St}_{X^\top X}$ can be written as

$$ \xi = X(X^\top X)^{-1}\Omega + \Delta \text{ with } \Omega \in \text{Skew}(p), $$
any matrix $Z \in \mathbb{R}^{n,p}$ can be naturally decomposed as

$$
\begin{align*}
Z &= X(X^\top X)^{-1} X^\top Z + (I_n - X(X^\top X)^{-1} X^\top) Z \\
&= X(X^\top X)^{-1} \operatorname{skew}(X^\top Z) + (I_n - X(X^\top X)^{-1} X^\top) Z + X(X^\top X)^{-1} \operatorname{sym}(X^\top Z) \\
&= \operatorname{Proj}_X[Z] + \operatorname{Proj}_X^{\perp}[Z],
\end{align*}
$$

where $\operatorname{Proj}_X$ and $\operatorname{Proj}_X^{\perp}$ are the operators defined by

$$
(7.13) \quad \operatorname{Proj}_X[Z] := X(X^\top X)^{-1} \operatorname{skew}(X^\top Z) + (I_n - X(X^\top X)^{-1} X^\top) Z,
$$

$$
(7.14) \quad \operatorname{Proj}_X^{\perp}[Z] := \operatorname{Id}_{\mathcal{E}} - \operatorname{Proj}_X(Z) = X(X^\top X)^{-1} \operatorname{sym}(X^\top Z).
$$

It is straightforward to verify that $\operatorname{Proj}_X$ and $\operatorname{Proj}_X^{\perp}$ are linear (oblique) projectors. In view of (7.14), it
is clear that this choice of projectors corresponds to attaching to $T_X \operatorname{St}_{X^\top X}$ the normal space

$$
(7.15) \qquad \mathrm{N}_X \mathrm{St}_{X^\top X} := \mathrm{Range}(\mathrm{Proj}_X^{\perp}) = \{X(X^\top X)^{-1} S \mid S \in \mathrm{Sym}(p)\},
$$

to be compared with (7.8). Given any symmetric operators $G_T(X) : \mathcal{E} \to \mathcal{E}$ and $G_N(X) : \mathcal{E} \to \mathcal{E}$,
positive definite on respectively $T_X \mathrm{St}_{X^\top X}$ and $\mathrm{N}_X \mathrm{St}_{X^\top X}$, $\mathrm{Proj}_X$ and $\mathrm{Proj}_X^\perp$ are orthogonal projectors
for the metric (7.12) with

$$
\widetilde{G}_T(X) = \mathrm{Proj}_X^{*,\mathcal{E}} G_T(X) \mathrm{Proj}_X, \quad \widetilde{G}_N(X) = (\mathrm{Proj}_X^{\perp})^{*,\mathcal{E}} G_N(X) \mathrm{Proj}_X^{\perp}.
$$

The following proposition provides expressions for the Euclidean adjoints of the tangent projection $\operatorname{Proj}_X$ and $\operatorname{Proj}_X^{\perp}$.

**Proposition 7.8.** The adjoint operators $\operatorname{Proj}_X^{*,\mathcal{E}}$ and $(\operatorname{Proj}_X^{\perp})^{*,\mathcal{E}}$ of the projection operators $\operatorname{Proj}_X$ and $\operatorname{Proj}_X^{\perp}$—with respect to the Euclidean inner product—are given by

$$
(7.16) \quad \operatorname{Proj}_X^{*,\mathcal{E}}[Z] = X \operatorname{skew}((X^\top X)^{-1} X^\top Z) + (I_n - X(X^\top X)^{-1} X^\top) Z, \quad Z \in \mathbb{R}^{n \times p},
$$

$$
(7.17) \qquad (\operatorname{Proj}_X^{\perp})^{*,\mathcal{E}}[Z] = X \operatorname{sym}((X^\top X)^{-1} X^\top Z), \quad Z \in \mathbb{R}^{n \times p}.
$$

In particular, the Euclidean orthogonal spaces of the tangent and normal spaces are

$$
(7.18) \quad (T_X \mathrm{St}_{X^\top X})^{\perp, \mathcal{E}} = \mathrm{Range}((\mathrm{Proj}_X^{\perp})^{*, \mathcal{E}}) = \{XS \mid S \in \mathrm{Sym}(p)\},
$$

$$
(7.19) \quad (\mathrm{N}_X \mathrm{St}_{X^\top X})^{\perp, \mathcal{E}} = \mathrm{Range}(\mathrm{Proj}_X^{*,\mathcal{E}}) = \{X\Omega + \Delta \mid \Omega \in \mathrm{Skew}(p) \text{ and } \Delta \in \mathbb{R}^{n \times p} \text{ with } X^\top \Delta = 0\}.
$$

Proof. By using the self-adjointness of the operators $X(X^\top X)^{-1}X^\top$ and skew, we have for any $\xi, \zeta \in \mathbb{R}^{n \times p}$,

$$
\begin{align*}
\langle \operatorname{Proj}_X[\xi], \zeta \rangle &= \langle X(X^\top X)^{-1} \operatorname{skew}(X^\top \xi), \zeta \rangle + \langle (I_n - X(X^\top X)^{-1} X^\top) \xi, \zeta \rangle \\
&= \langle \operatorname{skew}(X^\top \xi), (X^\top X)^{-1} X^\top \zeta \rangle + \langle \xi, (I_n - X(X^\top X)^{-1} X^\top) \zeta \rangle \\
&= \langle X^\top \xi, \operatorname{skew}((X^\top X)^{-1} X^\top \zeta) \rangle + \langle \xi, (I_n - X(X^\top X)^{-1} X^\top) \zeta \rangle \\
&= \langle \xi, X \operatorname{skew}((X^\top X)^{-1} X^\top \zeta) \rangle + (I_n - X(X^\top X)^{-1} X^\top) \zeta.
\end{align*}
$$

Formula (7.17) can be obtained identically, or from $(\operatorname{Proj}_X^{\perp})^{*,\mathcal{E}} = \operatorname{Id}_{\mathcal{E}} - \operatorname{Proj}_X^{*,\mathcal{E}}$. $\square$

<figure><img src="image_2.png" /></figure>

We mention the following characterization of the adjoint of the operator $\mathrm{Dc}(X)$ with respect to the
metric $g$ defined in (7.12).

**Proposition 7.9.** For any $S \in \mathrm{Sym}(p)$,

$$
\mathrm{Dc}(X)^{*,g}[S] = \widetilde{G}_N(X)^{-1}[XS],
$$

where $\widetilde{G}_N(X)^{-1} : (T_X \mathrm{St}_{X^\top X})^{\perp, \mathcal{E}} \to \mathrm{N}_X \mathrm{St}_{X^\top X}$ is the inverse of the operator $\widetilde{G}_N$ defined in (4.8).

*Proof.* Due to (3.12) and proposition 7.2:

$$
(7.20) \qquad \mathrm{Dc}(X)^{*,g}[S] = G(X)^{-1}\mathrm{Dc}(X)^{*,\mathcal{E}}[S] = G(X)^{-1}[XS].
$$

Then, in view of proposition 7.8, $\operatorname{Proj}_X^{*,\mathcal{E}}[XS] = 0$ and $(\operatorname{Proj}_X^{\perp})^{*,\mathcal{E}}[XS] = XS$. The result follows by applying (4.9) to (7.20). $\square$

<figure><img src="image_3.png" /></figure>
From section 4.1, the tangent and normal steps of the landing algorithm (4.12) read

$$ (7.21) \qquad d_T(X) = -\tilde{G}_T(X)^{-1} \mathrm{Proj}_X^{*,\mathcal{E}}[\nabla_{\mathcal{E}} f(X)], $$

$$ (7.22) \qquad d_N(X) = -\mathrm{Proj}_X^{\perp} \mathrm{Dc}(X)^{\dagger, \mathcal{E}} H(X)[c(X)]. $$

The operators $\tilde{G}_T(X)^{-1}$ and $H(X)$ can be freely chosen by the user which gives some latitude for the computation of $d_T(X)$ and $d_N(X)$, and which have rather clear physical interpretations. We obtain the following formulas for $d_T(X)$ and $d_N(X)$.

**Proposition 7.10.** The tangent and normal terms (eq. (7.21) and (7.22)) of the landing algorithm (1.3) relatively to the metric (7.12) are given by

$$ (7.23) \qquad d_T(X) = -\tilde{G}_T(X)^{-1}[X \operatorname{skew}((X^\top X)^{-1} X^\top \nabla_{\mathcal{E}} f(X)) + (I_n - X(X^\top X)^{-1} X^\top) \nabla_{\mathcal{E}} f(X)], $$

$$ (7.24) \qquad d_N(X) = -X(X^\top X)^{-1} \operatorname{sym}(X^\top X S), $$

where $S$ is the solution to the Sylvester equation (7.9) with $T = \frac{1}{2}H(X)[X^\top X - I_p]$. The expression of $d_N(X)$ becomes explicit in the following cases:

(i) if $H(X) = \mathrm{Id}_{\mathcal{E}}$, it reads $d_N(X) = -\frac{1}{2}X(I_p - (X^\top X)^{-1})$;

(ii) if $H(X) = \mathrm{Dc}(X)\mathrm{Dc}(X)^{*,g}$, it reads $d_N(X) = -\frac{1}{2}\tilde{G}_N(X)^{-1}[X(X^\top X - I_p)]$;

(iii) if $H(X) = \mathrm{Dc}(X)\mathrm{Dc}(X)^{*,\mathcal{E}}$, it reads $d_N(X) = -X(X^\top X - I_p)$.

*Proof.* Formula (7.23) follows immediately from (7.21) and (7.16). Formula (7.24) is obtained by combining (7.22) and the result of (7.6). For the cases (i) and (iii), we observe that due to (7.22),

$$ d_N(X) = \mathrm{Proj}_X^{\perp}[d_N(X)^{\mathcal{E}}] $$

where $d_N(X)^{\mathcal{E}}$ is the corresponding normal direction obtained in proposition 7.7. Since

$$ X(I_p - (X^\top X)^{-1}) = X(X^\top X)^{-1}(X^\top X - I_p) \in N_X \mathrm{St}_{X^\top X}, $$

$$ X(X^\top X - I_p) = X(X^\top X)^{-1}((X^\top X)^2 - X^\top X) \in N_X \mathrm{St}_{X^\top X}, $$

we have in both cases $d_N(X) = d_N(X)^{\mathcal{E}}$, which proves (i) and (iii). Finally, $H(X) = \mathrm{Dc}(X)\mathrm{Dc}(X)^{*,g}$ implies that $d_N(X) = -\mathrm{Dc}(X)^{*,g}c(X)$, which yields (ii) after using proposition 7.9. $\square$

It remains to choose the operator $\tilde{G}_T(X)^{-1}$ in (7.23), and the operator $\tilde{G}_N(X)^{-1}$ if one chooses (ii) for the normal direction.

### 7.2.2. Construction of the tangent and normal steps based on canonical tangent and normal metric representers

There is a natural choice of operators $\tilde{G}_T(X) : T_X \mathrm{St}_{X^\top X} \to N_X \mathrm{St}_{X^\top X}^{\perp, \mathcal{E}}$ and $\tilde{G}_N(X) : N_X \mathrm{St}_{X^\top X} \to T_X \mathrm{St}_{X^\top X}^{\perp, \mathcal{E}}$ with explicit inverses $\tilde{G}_T(X)^{-1}$ and $\tilde{G}_N(X)^{-1}$, leading to explicit expressions for the tangent and normal steps (7.23) and (7.24). Recalling (7.4) and (7.19), and (7.15) and (7.18), these ‘canonical’ operators read

$$ (7.25) \qquad \tilde{G}_T(X) : \quad T_X \mathrm{St}_{X^\top X} \to N_X \mathrm{St}_{X^\top X}^{\perp, \mathcal{E}} \\ X(X^\top X)^{-1} \Omega + \Delta \mapsto X \Omega + \Delta, $$

$$ (7.26) \qquad \tilde{G}_N(X) : \quad N_X \mathrm{St}_{X^\top X} \to T_X \mathrm{St}_{X^\top X}^{\perp, \mathcal{E}} \\ X(X^\top X)^{-1} S \mapsto X S, $$

where $\Omega \in \mathrm{Skew}(p)$, $\Delta \in \mathbb{R}^{n \times p}$ with $\Delta^\top X = 0$ and $S \in \mathrm{Sym}(p)$. These mappings have the following natural extensions to the whole $\mathcal{E} = \mathbb{R}^{n \times p}$:

$$ (7.27) \qquad G_T(X)[Z] := XX^\top Z + (I_n - X(X^\top X)^{-1} X^\top)Z, \quad Z \in \mathbb{R}^{n \times p}, $$

$$ (7.28) \qquad G_N(X)[Z] := XX^\top Z, \quad Z \in \mathbb{R}^{n \times p}. $$

With these definitions, we have indeed $G_T(X)\mathrm{Proj}_X = \tilde{G}_T(X)\mathrm{Proj}_X$ and $G_N(X)\mathrm{Proj}_X = \tilde{G}_N(X)\mathrm{Proj}_X$.

**Proposition 7.11.** The operators $G_T(X)$ and $G_N(X)$ of (7.27) and (7.28) are symmetric positive operators on $\mathcal{E}$, definite respectively on $\mathcal{E}$ and $N_X \mathrm{St}_{X^\top X}$. Formula (7.12) defines thus the associated metric

$$ (7.29) \qquad \begin{aligned} g(\xi, \zeta) &= \langle \xi, \mathrm{Proj}_X^{*,\mathcal{E}} G_T(X) \mathrm{Proj}_X \zeta \rangle + \langle \xi, (\mathrm{Proj}_X^{\perp})^{*,\mathcal{E}} G_N(X) \mathrm{Proj}_X \zeta \rangle \\ &= \langle \xi, (XX^\top + I_n - X(X^\top X)^{-1} X^\top) \zeta \rangle. \end{aligned} $$
Proof. The symmetry, the positivity of $G_T$ and $G_N$ is clear, as well as the definiteness of $G_T$ on $\mathcal{E}$. The definiteness of $G_N$ on $N_X \text{St}_{X^\top X}$ follows from the inequality

$$ \langle X(X^\top X)^{-1} S, G_N(X) X(X^\top X)^{-1} S \rangle = \langle X(X^\top X)^{-1} S, X S \rangle = \langle S, S \rangle, \quad \forall S \in \text{Sym}(p). $$

We then find that

$$ \begin{aligned} \text{Proj}_X^{*,\mathcal{E}} G_T(X) \text{Proj}_X \zeta &= \text{Proj}_X^{*,\mathcal{E}} [X \text{ skew}(X^\top \zeta) + (I_n - X(X^\top X)^{-1} X^\top) \zeta] \\ &= X \text{ skew}(X^\top \zeta) + (I_n - X(X^\top X)^{-1} X^\top) \zeta, \\ (\text{Proj}_X^{\perp})^{*,\mathcal{E}} G_N(X) \text{Proj}_X \zeta &= (\text{Proj}_X^{\perp})^{*,\mathcal{E}} [X \text{ sym}(X^\top \zeta)] = X \text{ sym}(X^\top \zeta). \end{aligned} $$

With these formulas, we obtain thus

$$ \begin{aligned} g(\xi, \zeta) &= \langle \xi, X \text{ skew}(X^\top \zeta) + (I_n - X(X^\top X)^{-1} X^\top) \zeta \rangle + \langle \xi, X \text{ sym}(X^\top \zeta) \rangle \\ &= \langle \xi, (XX^\top + I_n - X(X^\top X)^{-1} X^\top) \zeta \rangle. \end{aligned} $$

□

Since $\widetilde{G}_T^{-1}$ and $\widetilde{G}_N^{-1}$ are explicit, we can provide explicit formulas for the tangent step $d_T(X) = -\text{grad}_{\text{St}_{X^\top X}} f(X)$ and the normal step $d_N(X) = -\nabla_g \psi(X)$ (corresponding to $H(X) = \text{Dc}(X)\text{Dc}(X)^{*,g}$) in the metric $g$ of (7.29).

**Proposition 7.12.** With the metric $g$ defined in (7.29), the tangent and normal directions of the landing algorithm (1.3) read

$$ (7.30) \qquad d_T(X) = -X(X^\top X)^{-1} \text{skew}((X^\top X)^{-1} X^\top \nabla_\mathcal{E} f(X)) - (I_n - X(X^\top X)^{-1} X^\top) \nabla_\mathcal{E} f(X), $$

$$ (7.31) \qquad d_N(X) = -\frac{1}{2} X(I_p - (X^\top X)^{-1}). $$

*Proof.* The inverse of the operators $\widetilde{G}_T(X)$ and $\widetilde{G}_N(X)$ read obviously

$$ (7.32) \qquad \widetilde{G}_T(X)^{-1} : \begin{array}{ccc} N_X \text{St}_{X^\top X}^{\perp, \mathcal{E}} & \longrightarrow & T_X \text{St}_{X^\top X} \\ X \Omega + \Delta & \longmapsto & X(X^\top X)^{-1} \Omega + \Delta, \end{array} $$

$$ (7.33) \qquad \widetilde{G}_N(X)^{-1} : \begin{array}{ccc} T_X \text{St}_{X^\top X}^{\perp, \mathcal{E}} & \longrightarrow & N_X \text{St}_{X^\top X} \\ X S & \longmapsto & X(X^\top X)^{-1} S. \end{array} $$

Formula (7.30) follow from (7.23), and combining (7.33) with the item (ii) of proposition 7.10 yields

$$ d_N(X) = -\frac{1}{2} X(X^\top X)^{-1} (X^\top X - I_p) = -\frac{1}{2} X(I_p - (X^\top X)^{-1}). $$

□

**Remark 7.1.** Here, $d_N(X)$ is both the Riemannian gradient of the penalty function $\psi(X)$ and the 'Newton-like' direction (7.22) with $H(X) = \text{Id}_\mathcal{E}$. This happens, according to proposition 4.4, when $\widetilde{G}_N(X) = \text{Dc}(X)^{*,\mathcal{E}} \text{Dc}(X)$ as operators from $N_X \text{St}_{X^\top X}$ to $T_X \text{St}_{X^\top X}^{\perp, \mathcal{E}}$. This is incidentally indeed the case here: for any $S \in \text{Sym}(p)$,

$$ \begin{aligned} \text{Dc}(X)^{*,\mathcal{E}} \text{Dc}(X)[X(X^\top X)^{-1} S] &= \text{Dc}(X)^{*,\mathcal{E}} [\text{sym}(X^\top X(X^\top X)^{-1} S)] = \text{Dc}(X)^{*,\mathcal{E}} [S] \\ &= XS = \widetilde{G}_N(X)[X(X^\top X)^{-1} S]. \end{aligned} $$

### 7.2.3. Construction of the tangent and normal steps based on the $\beta$-metric

It turns out that the projectors $\text{Proj}_X$ and $\text{Proj}_X^{\perp}$ of (7.13) and (7.14) are exactly the orthogonal projectors for the family of $\beta$-metric $g^\beta$ considered in Goyens et al. (2026):

$$ (7.34) \qquad g^\beta(\xi, \zeta) = \langle (I - (1 - \beta)X(X^\top X)^{-1} X^\top) \zeta(X^\top X)^{-1}, \xi \rangle, $$

see Goyens et al. (2026), Proposition 4. The $\beta$-metric $g^\beta$ is the pull-back of a generalization of the canonical metric Edelman et al. (1998) on the Stiefel manifold $\text{St}$ to the layered manifold $\text{St}_{X^\top X}$. We observe, however, that $g^\beta(\xi, \zeta)$ does not coincide with the metric $g$ found in (7.29), despite both generate the same family of normal spaces $N_X \text{St}_{X^\top X}$.

In the next proposition, we show that the $\beta$-metric (7.34) is a special case of (7.12) for particular choices of $G_T(X)$ and $G_N(X)$. Then, we verify that the operators $\widetilde{G}_T$ and $\widetilde{G}_N$ have explicit inverses, which is the reason why this family of metrics leads to explicit formulas for the tangent and normal terms $d_T(X)$ and $d_N(X)$.
**Proposition 7.13.** Let $\beta > 0$. The $\beta$-metric (7.34) is the metric (7.12) with $G_T(X) : \mathcal{E} \to \mathcal{E}$ and $G_N(X) : \mathcal{E} \to \mathcal{E}$ being the symmetric operators (with respect to the Euclidean inner product) defined by

$$ (7.35) \qquad G_T(X)[\xi] := ((I_n - X(X^\top X)^{-1}X^\top)\xi + \beta X(X^\top X)^{-1}X^\top\xi)(X^\top X)^{-1}, $$

$$ (7.36) \qquad G_N(X)[\xi] := \beta X(X^\top X)^{-1}X^\top\xi(X^\top X)^{-1}. $$

The operator $G_T(X)$ is positive definite on $\mathbb{R}^{n \times p}$ and the operator $G_N(X)$ is positive definite on $N_X \text{St}_{X^\top X}$.

*Proof.* Let us denote by $\Pi_X := X(X^\top X)^{-1}X^\top$ the projection matrix onto the image space of $X$. We observe that for $\xi, \zeta \in \mathbb{R}^{n \times p}$,

$$ g^{\beta}(\xi, \zeta) = \beta \langle \Pi_X \xi (X^\top X)^{-1}, \Pi_X \zeta \rangle + \langle (I_n - \Pi_X) \xi (X^\top X)^{-1}, (I_n - \Pi_X) \zeta \rangle. $$

Consequently,

$$ \begin{aligned} g^{\beta}(\text{Proj}_X[\xi], \text{Proj}_X^{\perp}[\zeta]) &= \langle X(X^\top X)^{-1} \text{skew}(X^\top \xi)(X^\top X)^{-1}, X(X^\top X)^{-1} \text{sym}(X^\top \zeta) \rangle \\ &= \langle \text{skew}(X^\top \xi), (X^\top X)^{-1} \text{sym}(X^\top \zeta)(X^\top X)^{-1} \rangle = 0, \end{aligned} $$

due to the Frobenius orthogonality between skew and symmetric matrices. We have thus retrieved the fact that the tangent space $T_X \text{St}_{X^\top X}$ and the normal space $N_X \text{St}_{X^\top X}$ of (7.15) are orthogonal for the metric $g^\beta$. This implies in particular that

$$ g^{\beta}(\xi, \zeta) = g^{\beta}(\text{Proj}_X[\xi], \text{Proj}_X[\zeta]) + g^{\beta}(\text{Proj}_X^{\perp}[\xi], \text{Proj}_X^{\perp}[\zeta]), $$

with

$$ \begin{aligned} g^{\beta}(\text{Proj}_X[\xi], \text{Proj}_X[\zeta]) &= \beta \langle \Pi_X \text{Proj}_X[\xi](X^\top X)^{-1}, \text{Proj}_X[\zeta] \rangle + \langle (I_n - \Pi_X) \text{Proj}_X[\xi](X^\top X)^{-1}, \text{Proj}_X[\zeta] \rangle \\ &= \langle G_T(X) \text{Proj}_X[\xi], \text{Proj}_X[\zeta] \rangle, \end{aligned} $$

$$ \begin{aligned} g^{\beta}(\text{Proj}_X^{\perp}[\xi], \text{Proj}_X^{\perp}[\zeta]) &= \beta \langle \Pi_X \text{Proj}_X[\xi](X^\top X)^{-1}, \text{Proj}_X[\zeta] \rangle_{\mathcal{E}} \\ &= \langle G_N(X) \text{Proj}_X[\xi], \text{Proj}_X[\zeta] \rangle. \end{aligned} $$

Thus, $g^\beta$ is the metric $g$ of (7.12) for $G_T$ and $G_N$ defined by (7.35) and (7.36). The positive definiteness of $G_T(X)$ and $G_N(X)$ are visible from

$$ \begin{aligned} \langle G_T(X)[\xi], \xi \rangle &= \beta \|\Pi_X \xi (X^\top X)^{-\frac{1}{2}}\|^2 + \|\langle (I_n - \Pi_X) \xi (X^\top X)^{-\frac{1}{2}} \rangle\|^2, \\ \forall S \in \text{Sym}(p), \quad \langle G_N(X)[X(X^\top X)^{-1}S], X(X^\top X)^{-1}S \rangle \\ &= \beta \langle X(X^\top X)^{-1}S(X^\top X)^{-1}, X(X^\top X)^{-1}S \rangle \\ &= \beta \|(X^\top X)^{-\frac{1}{2}} S(X^\top X)^{-\frac{1}{2}}\|^2. \end{aligned} $$

□

**Proposition 7.14.** The inverse of the operators

$$ \tilde{G}_T(X) = \text{Proj}_X^{*,\mathcal{E}} G_T(X) \text{Proj}_X \text{ and } \tilde{G}_N(X) = (\text{Proj}_X^{\perp})^{*,\mathcal{E}} G_N(X) \text{Proj}_X^{\perp}, $$

associated to the $\beta$-metric $g^\beta$ are explicit and given by

$$ (7.37) \qquad \tilde{G}_T(X)^{-1} : \begin{array}{ccc} N_X \text{St}_{X^\top X}^{\perp, \mathcal{E}} & \longrightarrow & T_X \text{St}_{X^\top X} \\ X \Omega + \Delta & \longmapsto & \frac{1}{\beta} X \Omega X^\top X + \Delta X^\top X, \end{array} $$

$$ (7.38) \qquad \tilde{G}_N(X)^{-1} : \begin{array}{ccc} T_X \text{St}_{X^\top X}^{\perp, \mathcal{E}} & \longrightarrow & N_X \text{St}_{X^\top X} \\ X S & \longmapsto & \frac{1}{\beta} X S(X^\top X), \end{array} $$

for any $\Omega \in \text{Skew}(p)$, $\Delta \in \mathbb{R}^{n \times p}$ with $\Delta^\top X = 0$ and $S \in \text{Sym}(p)$. Consequently, the tangent and normal directions of the landing algorithm (1.3) with $H(X) = \text{Dc}(X)\text{Dc}(X)^{*;g^\beta}$ read

$$ \begin{aligned} d_T(X) &= -\text{grad}_{\text{St}_{X^\top X}}^{g^\beta} f(X) \\ (7.39) \qquad &= -\frac{1}{\beta} X \text{ skew}((X^\top X)^{-1}X^\top \nabla_{\mathcal{E}} f(X)) X^\top X - (I_n - X(X^\top X)^{-1}X^\top) \nabla_{\mathcal{E}} f(X) X^\top X, \end{aligned} $$

$$ (7.40) \qquad d_N(X) = -\nabla_{g^\beta} \psi(X) = -\frac{1}{2\beta} X(X^\top X - I_p)X^\top X. $$
Proof. Using (7.35) and (7.36), we see that $\tilde{G}_T(X)$ and $\tilde{G}_N(X)$ are the operators

$$ \tilde{G}_T(X) : \quad T_X \mathrm{St}_{X^\top X} \quad \rightarrow \quad N_X \mathrm{St}_{X^\top X}^{\perp, \varepsilon} \\ X(X^\top X)^{-1}\Omega + \Delta \quad \longmapsto \quad \beta X(X^\top X)^{-1}\Omega(X^\top X)^{-1} + \Delta(X^\top X)^{-1}, $$

$$ \tilde{G}_N(X) : \quad N_X \mathrm{St}_{X^\top X} \quad \rightarrow \quad T_X \mathrm{St}_{X^\top X}^{\perp, \varepsilon} \\ X(X^\top X)^{-1} S \quad \longmapsto \quad \beta X(X^\top X)^{-1} S(X^\top X)^{-1}, $$

whose inverses are obviously given by (7.37) and (7.38). Inserting these formulas into (7.21) and item (ii) of proposition 7.10 yields the result. □

The reader may verify that (7.39) coincide with the formula derived in Goyens et al. (2026). In this latter work, the normal term was considered to be the Euclidean gradient $d_N(X) = -\nabla_{\varepsilon}\psi(X)$, which corresponds to (7.22) with $H(X) = Dc(X)Dc(X)^{*,*}$ (item (iii) of proposition 7.10).

REFERENCES

Ablin, P. and Peyré, G. (2022). Fast and accurate optimization on the orthogonal manifold without retraction. In *International Conference on Artificial Intelligence and Statistics*, pages 5636–5657. PMLR.

Absil, P.-A., Mahony, R., and Sepulchre, R. (2008). *Optimization Algorithms on Matrix Manifolds*. Princeton University Press.

Absil, P.-A. and Malick, J. (2012). Projection-like Retractions on Matrix Manifolds. *SIAM Journal on Optimization*, 22(1):135–158.

Absil, P.-A., Trumpf, J., Mahony, R., and Andrews, B. (2009). All roads lead to Newton: Feasible second-order methods for equality-constrained optimization. https://sites.uclouvain.be/absil/2009.024.

Allaire, G., Dapogny, C., and Jouve, F. (2021). Chapter 1 - Shape and topology optimization. In Bonito, A. and Nochetto, R. H., editors, *Handbook of Numerical Analysis*, volume 22 of *Geometric Partial Differential Equations - Part II*, pages 1–132. Elsevier.

Bai, Y. and Mei, S. (2018). Analysis of sequential quadratic programming through the lens of riemannian optimization. *arXiv preprint arXiv:1805.08756*.

Barbarosie, C., Toader, A.-M., and Lopes, S. (2020). A gradient-type algorithm for constrained optimization with application to microstructure optimization. *Discrete and Continuous Dynamical Systems - B*, 25(5):1729–1755.

Berahas, A. S., Curtis, F. E., O'Neill, M. J., and Robinson, D. P. (2024). A Stochastic Sequential Quadratic Optimization Algorithm for Nonlinear-Equality-Constrained Optimization with Rank-Deficient Jacobians. *Mathematics of Operations Research*, 49(4):2212–2248.

Berahas, A. S., Curtis, F. E., Robinson, D., and Zhou, B. (2021). Sequential quadratic optimization for nonlinear equality constrained stochastic optimization. *SIAM Journal on Optimization*, 31(2):1352–1379.

Biros, G. and Ghattas, O. (2005). Parallel Lagrange–Newton–Krylov–Schur Methods for PDE-Constrained Optimization. Part I: The Krylov–Schur Solver. *SIAM Journal on Scientific Computing*, 27(2):687–713.

Boggs, P. T. and Tolle, J. W. (1995). Sequential Quadratic Programming. *Acta Numerica*, 4:1–51.

Booker, A. and Ong, C.-y. (1971). Multiple-constraint adaptive filtering. *Geophysics*, 36(3):498–509.

Conn, A. R., Gould, N. I. M., and Toint, P. (1991). A Globally Convergent Augmented Lagrangian Algorithm for Optimization with General Constraints and Simple Bounds. *SIAM Journal on Numerical Analysis*, 28(2):545–572.

Curtis, F. E., O'Neill, M. J., and Robinson, D. P. (2024). Worst-case complexity of an sqp method for nonlinear equality constrained stochastic optimization. *Mathematical Programming*, 205(1):431–483.

Dietrich, A., Ott, C., and Albu-Schäffer, A. (2015). An overview of null space projections for redundant, torque-controlled robots. *The International Journal of Robotics Research*, 34(11):1385–1400.

Edelman, A., Arias, T. A., and Smith, S. T. (1998). The geometry of algorithms with orthogonality constraints. *SIAM journal on Matrix Analysis and Applications*, 20(2):303–353.

Evtushenko, Y. G. and Zhadan, V. G. (1994). Stable barrier-projection and barrier-Newton methods in linear programming. *Computational Optimization and Applications*, 3(4):289–303.

Fang, Y., Na, S., Mahoney, M. W., and Kolar, M. (2024). Fully Stochastic Trust-Region Sequential Quadratic Programming for Equality-Constrained Optimization Problems. *SIAM Journal on Optimization*, 34(2):2007–2037.
Feppon, F. (2024). Density-based topology optimization with the Null Space Optimizer: A tutorial and a comparison. *Structural and Multidisciplinary Optimization*, 67(1).

Feppon, F., Allaire, G., and Dapogny, C. (2020a). Null space gradient flows for constrained optimization with applications to shape optimization. *ESAIM: Control, Optimisation and Calculus of Variations*, page 90.

Feppon, F., Allaire, G., Dapogny, C., and Jolivet, P. (2020b). Topology optimization of thermal fluid-structure systems using body-fitted meshes and parallel computing. *Journal of Computational Physics*, 417:109574.

Feppon, F. and Lermusiaux, P. F. J. (2019). The Extrinsic Geometry of Dynamical Systems Tracking Nonlinear Matrix Projections. *SIAM Journal on Matrix Analysis and Applications*, 40(2):814–844.

Fletcher, R., Gould, N. I M., Leyffer, S., Toint, P. L., and Wächter, A. (2002). Global Convergence of a Trust-Region SQP-Filter Algorithm for General Nonlinear Programming. *SIAM Journal on Optimization*, 13(3):635–659.

Frost, O. (1972). An algorithm for linearly constrained adaptive array processing. *Proceedings of the IEEE*, 60(8):926–935.

Gangi, A. F. and Byun, B. S. (1976). The corrective gradient projection method and related algorithms applied to seismic array processing. *Geophysics*, 41(5):970–984.

Gao, B., Vary, S., Ablin, P., and Absil, P. A. (2022). Optimization flows landing on the Stiefel manifold. *IFAC-PapersOnLine*, 55(30):25–30.

Gharaei, A., Amjadian, A., Amjadian, A., Shavandi, A., Hashemi, A., Taher, M., and Mohamadi, N. (2023). An integrated lot-sizing policy for the inventory management of constrained multi-level supply chains: Null-space method. *International Journal of Systems Science: Operations & Logistics*, 10(1):2083254.

Gill, P., Wong, E., Murray, W., and Saunders, M. (2015). User's Guide for SNOPT Version 7.4: Software for Large-Scale Nonlinear Programming.

Gill, P. E., Murray, W., and Saunders, M. A. (2002). SNOPT: An SQP Algorithm for Large-Scale Constrained Optimization. *SIAM Journal on Optimization*, 12(4):979–1006.

Gill, P. E. and Wong, E. (2012). Sequential Quadratic Programming Methods. In Lee, J. and Leyffer, S., editors, *Mixed Integer Nonlinear Programming*, pages 147–224, New York, NY. Springer.

Goyens, F., Absil, P.-A., and Feppon, F. (2026). Geometric Design of the Tangent Term in Landing Algorithms for Orthogonality Constraints. In Nielsen, F. and Barbaresco, F., editors, *Geometric Science of Information*, pages 133–141, Cham. Springer Nature Switzerland.

Goyens, F., Eftekhari, A., and Boumal, N. (2024). Computing Second-Order Points Under Equality Constraints: Revisiting Fletcher's Augmented Lagrangian. *Journal of Optimization Theory and Applications*, 201(3):1198–1228.

Henrot, A. and Pierre, M. (2018). *Shape Variation and Optimization*, volume 28 of *EMS Tracts in Mathematics*. European Mathematical Society (EMS), Zürich.

Hu, X., Cheng, Y., Yang, D., Xu, Z., Yuan, Z., Yu, J., Xu, C., Jiang, Z., and Zhou, S. (2025). Ostquant: Refining large language model quantization with orthogonal and scaling transformations for better distribution fitting. *arXiv preprint arXiv:2501.13987*.

Hüper, K., Markina, I., and Leite, F. S. (2021). A Lagrangian approach to extremal curves on Stiefel manifolds. *Journal of Geometric Mechanics (JGM)*.

Jongen, H. T. and Stein, O. (2003). On the Complexity of Equalizing Inequalities. *Journal of Global Optimization*, 27(4):367–374.

Jongen, H. Th. and Stein, O. (2004). Constrained Global Optimization: Adaptive Gradient Flows. In Floudas, C. A. and Pardalos, P., editors, *Frontiers in Global Optimization*, pages 223–236, Boston, MA. Springer US.

Li, Z.-Y. and Zhou, B. (2018). Spectral decomposition based solutions to the matrix equation. *IET Control Theory & Applications*, 12(1):119–128.

Liu, X. and Yuan, Y. (2010). A null-space primal-dual interior-point algorithm for nonlinear optimization with nice convergence properties. *Mathematical programming*, 125(1):163–193.

Luenberger, D. G. (1972). The Gradient Projection Method Along Geodesics. *Management Science*, 18(11):620–631.

Mataigne, S., Zimmermann, R., and Miolane, N. (2025). An Efficient Algorithm for the Riemannian Logarithm on the Stiefel Manifold for a Family of Riemannian Metrics. *SIAM Journal on Matrix Analysis and Applications*, 46(2):879–905.
Miller, S. A. and Malick, J. (2005). Newton methods for nonsmooth convex minimization: Connections among -Lagrangian, Riemannian Newton and SQP methods. *Mathematical Programming*, 104(2):609-633.

Mishra, B. and Sepulchre, R. (2016). Riemannian preconditioning. *SIAM Journal on Optimization*, 26(1):635-660.

Nie, P.-y. (2004). A null space method for solving system of equations. *Applied Mathematics and computation*, 149(1):215-226.

Nocedal, J. and Overton, M. L. (1985). Projected Hessian Updating Algorithms for Nonlinearly Constrained Optimization. *SIAM Journal on Numerical Analysis*, 22(5):821-850.

Nocedal, J., Öztoprak, F., and Waltz, R. A. (2014). An interior point method for nonlinear programming with infeasibility detection capabilities. *Optimization Methods and Software*, 29(4):837-854.

Nocedal, J. and Wright, S. J. (2006). *Numerical Optimization*. Springer.

Obara, M., Okuno, T., and Takeda, A. (2022). Sequential Quadratic Optimization for Nonlinear Optimization Problems on Riemannian Manifolds. *SIAM Journal on Optimization*, 32(2):822-853.

Rosen, J. B. (1960). The Gradient Projection Method for Nonlinear Programming. Part I. Linear Constraints. *Journal of the Society for Industrial and Applied Mathematics*, 8(1):181-217.

Rosen, J. B. (1961). The Gradient Projection Method for Nonlinear Programming. Part II. Nonlinear Constraints. *Journal of the Society for Industrial and Applied Mathematics*, 9(4):514-532.

Schechtman, S., Tiapkin, D., Muehlebach, M., and Moulines, É. (2023). Orthogonal Directions Constrained Gradient Method: From non-linear equality constraints to Stiefel manifold. In *Proceedings of Thirty Sixth Conference on Learning Theory*, pages 1228-1258. PMLR.

Schropp, J. and Singer, I. (2000). A dynamical systems approach to constrained minimization. *Numerical Functional Analysis and Optimization*, 21(3-4):537-551.

Shi, Q. and Wang, X. (2025). Adaptive directional decomposition methods for nonconvex constrained optimization.

Shikhman, V. and Stein, O. (2009). Constrained Optimization: Projected Gradient Flows. *Journal of Optimization Theory and Applications*, 140(1):117-130.

Song, Y., Li, P., Gao, B., and Yuan, K. (2025). Distributed retraction-free and communication-efficient optimization on the stiefel manifold. *arXiv preprint arXiv:2506.02879*.

Tanabe, K. (1974). An algorithm for constrained maximization in nonlinear programming. *Journal of the Operations Research Society of Japan*, 17:184-201.

Tanabe, K. (1980). A geometric method in nonlinear programming. *Journal of Optimization Theory and Applications*, 30(2):181-210.

Tapia, R. A. (1977). Diagonalized multiplier methods and quasi-Newton methods for constrained optimization. *Journal of Optimization Theory and Applications*, 22(2):135-194.

Vary, S., Ablin, P., Gao, B., and Absil, P.-A. (2024). Optimization without Retraction on the Random Generalized Stiefel Manifold. In *Proceedings of the 41st International Conference on Machine Learning*, pages 49226-49248. PMLR.

Wächter, A. and Biegler, L. T. (2005). Line Search Filter Methods for Nonlinear Programming: Motivation and Global Convergence. *SIAM Journal on Optimization*, 16(1):1-31.

Wegert, Z. J., Roberts, A. P., and Challis, V. J. (2023). A Hilbertian projection method for constrained level set-based topology optimisation. *Structural and Multidisciplinary Optimization*, 66(9):204.

Yamashita, H. (1980). A differential equation approach to nonlinear programming. *Mathematical Programming*, 18(1):155-168.

Yuan, Y.-x. (2001). A null space algorithm for constrained optimization. *Advances in Scientific Computing, Science Press, Beijing*, pages 210-218.

Yulin, M. and Xiaoming, W. (2004). A level set method for structural topology optimization with multi-constraints and multi-materials. *Acta Mechanica Sinica*, 20(5):507-518.

Zhang, Y., Hu, J., Cui, J., Lin, L., Wen, Z., and Li, Q. (2024). Retraction-free optimization over the Stiefel manifold with application to the LoRA fine-tuning.