Department of Policy and Planning Sciences
Discussion Paper Series

No.1381

Riemannian Interior Point Methods
for Constrained Optimization on Manifolds
byZhijian LAI and Akiko YOSHISE

March 18, 2022
Revised December 4, 2022
Revised February 7, 2024

UNIVERSITY OF TSUKUBA
Tsukuba, Ibaraki 305-8573
JAPANRiemannian Interior Point Methods
for Constrained Optimization on Manifolds∗
Zhijian Lai†

Akiko Yoshise‡

March 18, 2022
Revised December 4, 2022
Revised February 7, 2024

Abstract
We extend the classical primal-dual interior point method from the Euclidean setting to the
Riemannian one. Our method, named the Riemannian interior point method, is for solving Riemannian constrained optimization problems. We establish its local superlinear and quadratic
convergence under the standard assumptions. Moreover, we show its global convergence when
it is combined with a classical line search. Our method is a generalization of the classical framework of primal-dual interior point methods for nonlinear nonconvex programming. Numericalexperiments show the stability and eﬃciency of our method.
Communicated by Sándor Zoltán Németh
Keywords: Riemannian manifolds, Riemannian optimization, Nonlinear optimization andInterior point method.

1

Introduction

In this paper, we consider the following Riemannian Constrained Optimization Problem (RCOP):
min

f pxq

s.t.

hpxq “ 0, and gpxq ď 0,

xPM

(RCOP)

where M is a connected, complete d-dimensional Riemannian manifold and f : M Ñ R, h : M Ñ Rl ,and g : M Ñ Rm are smooth functions. This problem appears in many applications, for instance,
matrix approximation with nonnegative constraints on a ﬁxed-rank manifold [31] and orthogonal
nonnegative matrix factorization on the Stiefel manifold [19]; for more, see [23, 27].
The body of knowledge on (RCOP) without h, g, often called simply Riemannian optimization,
has grown considerably in the last 20 years. The well-known methods in the Euclidean setting,such as steepest descent, Newton, and trust region, have been extended to the Riemannian setting[1, 16, 5]. By contrast, research on (RCOP) is still in its infancy. The earliest studies go backto ones on the optimality conditions. Yang et al. [36] extended the Karush Kuhn Tucker (KKT)
conditions to (RCOP). Bergmann and Herzog [2] considered more Constraint Qualiﬁcations (CQs)
on manifolds. Yamakawa and Sato [34] proposed sequential optimality conditions in the Riemanniancase. Liu and Boumal [23] were the ﬁrst to develop practical algorithms. They extended theaugmented Lagrangian method and exact penalty method to (RCOP). Schiela and Ortiz [30] andObara et al. [27] proposed the Riemannian sequential quadratic programming method. However,to our knowledge, interior point methods have yet to be considered for (RCOP).
The advent of interior point methods in the 1980s greatly advanced the ﬁeld of optimization[33, 37]. By the early 1990s, the success of these methods in linear and quadratic programming
ignited interest in using them on nonlinear nonconvex cases [11, 35]. From the 1990s to the ﬁrst
decade of the 21st century, a large number of interior point methods for nonlinear programming∗ An earlier version of this article has been circulated under the title “Superlinear and Quadratic Convergence ofRiemannian Interior Point Methods.”† University of Tsukuba, Tsukuba, Ibaraki 305-8573, Japan. galvin.lai@outlook.com‡ Corresponding author. University of Tsukuba, Tsukuba, Ibaraki 305-8573, Japan. yoshise@sk.tsukuba.ac.jp

1emerged. They proved to be as successful as the linear ones. A subclass known as primal-dual
interior point methods is the most eﬃcient and practical. As described in [24], the primal-dual
approach to linear programming was introduced in [25]: it was ﬁrst developed as an algorithm in[20] and eventually became standard for the nonlinear case as well [11, 35].
Contribution. In this paper, we extend the primal-dual interior point algorithms from the
Euclidean setting, i.e., M “ Rd in (RCOP), to the Riemannian setting. We call this extension the
Riemannian Interior Point Method (RIPM). Our contributions are summarized as follows:
1. To our knowledge, this is the ﬁrst study to apply the primal-dual interior point method tothe nonconvex constrained optimization problem on Riemannian manifolds. One signiﬁcant
contribution is that we establish many essential foundational concepts for the general interior
point method in the Riemannian context, such as the KKT vector ﬁeld and its covariant
derivative. In addition, we build the ﬁrst framework for the Riemannian version of the interior
point method. These contributions will have uses in the future, especially in developing moreadvanced interior point methods.
2. We give a detailed theoretical analysis to ensure local and global convergence of RIPM. Considering that many practical problems involve minimizing a nonconvex function on Riemannian
manifolds, the theoretical counterparts of our method are the early interior point methods fornonlinear nonconvex programming ﬁrst proposed by El-Bakry et al. [11].
3. Our numerical experiments1 demonstrate the great potential of RIPM. The method meets thechallenges presented in these experiments with better stability and higher accuracy comparedwith the current Riemannian exact penalty, augmented Lagrangian, and sequential quadratic
programming methods.
Organization. The rest of this paper is organized as follows. In Section 2, we review thenotation of Riemannian geometry and explain the Riemannian Newton method. In Section 3, we givea full interpretation of our RIPM and describe a prototype algorithm of RIPM. We also investigate
the use of Krylov subspace methods to eﬃciently solve a condensed form of a perturbed Newton
equation. This is particularly important for the numerical implementation of RIPM. Section 4 gives
the necessary preliminaries and auxiliary results needed to prove convergence in our subsequent
sections. In Section 5, we give the proof of local superlinear and quadratic convergence of the
prototype algorithm of RIPM. Section 6 describes a globally convergent version of RIPM with aclassical line search; then, Section 7 proves its global convergence. Section 8 is a collection ofnumerical experiments. Section 9 summarizes our research and presents future work.

2

Notation

2.1

Riemannian Geometry

Let us brieﬂy review some concepts from Riemannian geometry, following the notation of [5]. M
denotes a ﬁnite-dimensional smooth manifold. Let p P M and Tp M be the tangent space at p with
0p being its zero element. We use a canonical identiﬁcation Tp E – E for aŤvector space E and p P E.
A vector ﬁeld is a map V : M Ñ T M with V ppq P Tp M, where T M :“ pPM Tp M is the tangent
bundle. XpMq denotes the set of all smooth vector ﬁelds deﬁned on M. Furthermore, M is a
Riemannian manifold if it is equipped with a Riemannian metric, that is, a choice of inner productx¨, ¨yp : Tp M ˆ Tp M Ñ R for each tangent space at p on M such that for all V, W P XpMq, the mapp ÞÑ xV ppq, W ppqyp(1)a
is a smooth function from M to R. Riemannian metric induces the norm }ξ}p :“ xξ, ξyp forξ P Tp M. We often omit the subscript p if it is clear from the context. Throughout this paper, we
assume that all the manifolds involved are connected and complete. Given a curve segment on M,şb9
9 P Tcptq M is the velocity
c : ra, bs Ñ M, the length of c is deﬁned as Lpcq :“ a }cptq}
cptq dt, where cptq
vector of c at t. Since M is connected, there exists a curve segment connecting any pair of points
p, q P M. Indeed, M is a metric space under the Riemannian distance dpp, qq :“ inf c Lpcq where theinﬁmum is taken over all curve segments joining p to q [22, Thm. 2.55]. For two manifolds M1 , M2
1 The code is freely available at https://doi.org/10.5281/zenodo.10612799

2and a smooth map F : M1 Ñ M2 , the diﬀerential of F at p P M1 is a linear operator denoted
as DF ppq : Tp M1 Ñ TF ppq M2 . Let FpMq be the set of all smooth scalar ﬁelds (or say real-valued
functions) f : M Ñ R. The Riemannian gradient of f at p, grad f ppq, is deﬁned as the unique elementof Tp M that satisﬁes xξ, grad f ppqyp “ Df ppqrξs for all ξ P Tp M, where Df ppq : Tp M Ñ Tf ppq R – R.
Note that for any f P FpMq, the gradient vector ﬁeld x ÞÑ grad f pxq is a smooth vector ﬁeld on
M, i.e., grad f P XpMq. A retraction R : T M Ñ M is a smooth map such that Rp p0p q “ p and
DRp p0p q “ idTp M , i.e., the identity map on Tp M, where Rp is the restriction of R to Tp M andDRp p0p q : T0p pTp Mq – Tp M Ñ Tp M is the diﬀerential of Rp at 0p . One theoretically perfect typeof retraction is the exponential map, denoted as Exp. Since M is complete, the exponential map
is well-deﬁned on the whole tangent bundle. Let Expp : Tp M Ñ M be the exponential map at p;then t ÞÑ Expp ptξq is the unique geodesic that passes through p with velocity ξ P Tp M when t “ 0.

2.2

Riemannian Newton Method

The Newton method is a powerful tool for ﬁnding the zeros of nonlinear functions in the Euclidean
setting. The generalized Newton method has been studied in the Riemannian setting; it aims toﬁnd a singularity for the vector ﬁeld F P XpMq, i.e., a point p P M such that,F ppq “ 0p P Tp M.

(2)

Let ∇ be the Levi-Civita connection on M, i.e., the unique symmetric connection compatible with
the Riemannian metric. The covariant derivative ∇F assigns each point p P M a linear operator∇F ppq : Tp M Ñ Tp M. In particular, the Riemannian Hessian of f P FpMq at p is a self-adjoint
operator on Tp M, deﬁned as Hess f ppq :“ ∇pgrad f qppq. The Riemannian Newton method for solving(2) is performed as Algorithm 1.
Algorithm 1: Riemannian Newton Method for (2)Input: A vector ﬁeld F P XpMq, an initial point p0 P M and a retraction R on M.
Output: Sequence tpk u Ă M such that tpk u Ñ p˚ and F pp˚ q “ 0p˚ .
Set k Ñ 0;
while Stopping criterion not satisﬁed do
1. Obtain ξk P Tpk M by solving the Newton equation (a linear operator equation ontangent space Tpk M):∇F ppk qrξk s “ ´F ppk q;(3)
2. Compute the next point as pk`1 :“ Rpk pξk q;
3. k Ñ k ` 1;
end

3

Description of Riemannian Interior Point Methods

In this section, we will give a comprehensive interpretation of the Riemannian interior point method.
Following common usage in the interior-point literature, big letters denote the associated diagonal
matrix, e.g., Z “ diagpz1 , . . . , zm q with z P Rm . e denotes the all-ones vector, and 0 stands for zero
vector/matrix with proper dimensions.

3.1

KKT Vector Field and Its Covariant Derivative

The Lagrangian function of (RCOP) is
Lpx, y, zq :“ f pxq `

l
ÿ
i“1

yi hi pxq `mÿzi gi pxq,

i“1where y P Rl and z P Rm are Lagrange multipliers corresponding to the equality and inequality
constraints, respectively. With respect to the variable x, Lp¨, y, zq is a real-valued function deﬁned

3

on M, and its Riemannian gradient isgradx Lpx, y, zq “ grad f pxq `lÿ

yi grad hi pxq `

i“1

m
ÿ

zi grad gi pxq,

i“1where tgrad hi pxquli“1 and tgrad gi pxqum
i“1 are the gradients of the component functions of h and g.
The Riemannian KKT conditions (e.g., see [23, Deﬁnition 2.3] or [36]) for (RCOP) are given bygradx Lpx, y, zq “ 0x ; hpxq “ 0, gpxq ď 0, z ě 0; Zgpxq “ 0.
The above conditions can be written in terms of slack variables Rm Q s :“ ´gpxq, as
¨
˛
¨
˛
0x
gradx Lpx, y, zq
˚ hpxq‹˚‹‹ “ 0w “ ˚ 0 ‹ P T w N ,F pwq :“ ˚
˝ gpxq ` s
‚
˝ 0 ‚
ZSe
0

(4)

(5)

and pz, sq ě 0, where w :“ px, y, z, sq P N :“ M ˆ Rl ˆ Rm ˆ Rm . Here, we have generated a vector
ﬁeld F in (5) on the product manifold N , i.e., for any w P N , F pwq P Tw N – Tx M ˆ Rl ˆ Rm ˆ Rm .
We will call F : N Ñ T N in (5) above the KKT vector ﬁeld for (RCOP). Note that for ξ “pξx , ξy , ξz , ξs q P Tw N ,R̄w pξq :“ pRx pξx q, y ` ξy , z ` ξz , s ` ξs q
(6)
is a well-deﬁned retraction on the product manifold N as long as R is a retraction on M.
Now, we aim to ﬁnd a singularity of the KKT vector ﬁeld F on N under some nonnegative
constraints pz, sq ě 0. If the Riemannian Newton method is to be applied to (5), we must formulate
the covariant derivative of F at an arbitrary w P N . Before that, we need some new symbols. Fixing
a point x P M, we can deﬁne two linear operators Hx : Rl Ñ Tx M and Gx : Rm Ñ Tx M by
Hx rvs :“lÿ

Gx rus :“vi grad hi pxq,

i“1

m
ÿ

ui grad gi pxq,

(7)

i“1

respectively. Then, for Hx , its adjoint operator Hx˚ : Tx M Ñ Rl is given byT

Hx˚ rξs “ pxgrad h1 pxq, ξyx , ¨ ¨ ¨ , xgrad hl pxq, ξyx q .
TAlso, Gx˚ : Tx M Ñ Rm with Gx˚ rξs “ pxgrad g1 pxq, ξyx , ¨ ¨ ¨ , xgrad gm pxq, ξyx q .
Now, using the solutions of Exercises 5.4 and 5.13 on covariant derivatives of vector ﬁelds on
product manifolds in monograph [5], the following results can be easily derived. Given KKT vector
ﬁeld F in (5) and w P N , the covariant derivative of F at w is the operator ∇F pwq : Tw N Ñ Tw N
given by¨
˛
Hessx Lpwqr∆xs ` Hx r∆ys ` Gx r∆zs
˚ Hx˚ r∆xs‹‹,
∇F pwqr∆ws “ ˚(8)‚˝ Gx˚ r∆xs ` ∆sZ∆s ` S∆zwhere ∆w “ p∆x, ∆y, ∆z, ∆sq P Tw N and Hessx Lpwq is the Riemannian Hessian of real-valued
function Lp¨, y, zq. Indeed, Hessx Lpwq : Tx M Ñ Tx M satisﬁes
Hessx Lpwq “ Hess f pxq `lÿyi Hess hi pxq `

i“1

m
ÿ

zi Hess gi pxq,

(9)

i“1

where tHess hi pxquli“1 and tHess gi pxqum
i“1 are Hessians of the component functions of h and g.

3.2

Implication of Standard Assumptions

Next, we will discuss an important result about the covariant derivative of the KKT vector ﬁeld givenin (8). Let Apxq :“ ti | gi pxq “ 0u denote the active set at x P M. A prior study on Riemannian
optimality conditions [36, 2] showed that the following assumptions are meaningful. We call themthe standard Riemannian assumptions for (RCOP). Note that the x˚ and w˚ in (A2)-(A4) all referto those in (A1).
4

(A1) Existence. There exists w˚ “ px˚ , y ˚ , z ˚ , s˚ q satisfying the Riemannian KKT conditions (4).
Here, we introduce the slack variables s˚ :“ ´gpx˚ q.
(A2) Linear Independence Constraint Qualiﬁcation (LICQ).
Ťl
The set tgrad hi px˚ qui“1 tgrad gi px˚ quiPApx˚ q is linearly independent in Tx˚ M.
(A3) Strict complementarity. pz ˚ qi ą 0 if gi px˚ q “ 0 for all i “ 1, ¨ ¨ ¨ , m.
(A4) Second-order suﬃciency. xHessx Lpw˚ qrξs, ξy ą 0 for all nonzero ξ P Tx˚ Msatisfying xξ, grad hi px˚ qy “ 0 for i “ 1, . . . , l, and xξ, grad gi px˚ qy “ 0 for i P Apx˚ q.
Recall the Riemannian Newton method discussed in Section 2.2. It can be shown that if p˚ isa solution of equation (2) and the covariant derivative ∇F pp˚ q is nonsingular, then Algorithm 1
has the local superlinear convergence [12] and local quadratic convergence [13] under certain mild
conditions on the map p ÞÑ ∇F ppq. Note that the requirement of nonsingularity for the covariant
derivative at the solution point is of primary importance. Therefore, the next theorem motivates us
to use the Riemannian Newton method to solve (5).
Theorem 3.1. If the standard Riemannian assumptions (A1)-(A4) hold at some point w˚ , thenthe operator ∇F pw˚ q in (8) is nonsingular.
Proof. This proof omits all the asterisks of the variables. Deﬁne E :“ t1, . . . , lu and I :“ t1, . . . , mu.
Take some w “ px, y, z, sq P N satisfying (A1)-(A4), then we have si “ ´gi pxq and zi si “ 0for all i P I. For short, let A :“ Apxq Ă I. Suppose that ∇F pwqr∆ws “ 0 for some ∆w “
p∆x, ∆y, ∆z, ∆sq P Tw N – Tx M ˆ Rl ˆ Rm ˆ Rm . ∆yi denotes the components of the vector ∆y,as do ∆zi , ∆si . To prove its nonsingularity, we will show that ∆w “ 0. Expanding the equation∇F pwqr∆ws “ 0 givesÿÿ
$
0 “ Hessx Lpwqr∆xs `∆yi grad hi pxq ` ∆zi grad gi pxq,
’
’
’
’
iPE
iPI
’
&
0 “ xgrad hi pxq, ∆xy , for all i P E,
(10)
’
’
’
`
0
“
xgrad
g
pxq,
∆xy
∆s
,
for
all
i
P
I,
i
i
’
’
%
0 “ zi ∆si ` si ∆zi , for all i P I.
Strict complementarity (A3) and the last equalities above imply that ∆si “ 0 for all i P A and∆zi “ 0 for i P IzA. Substituting those values into the system (10) reduces it to
$
ÿ
ÿ
’
0
“
Hess
Lpwqr∆xs
`
∆y
grad
h
pxq
`
∆zi grad gi pxq,
x
i
i
’
’
&
iPE
iPA
(11)
0 “ xgrad hi pxq, ∆xy , for all i P E,’’’%
0 “ xgrad gi pxq, ∆xy , for all i P A,
and ∆si “ ´ xgrad gi pxq, ∆xy for all i P IzA. It follows from system (11) thatÿ
ÿ
0 “ xHessx Lpwqr∆xs `∆yi grad hi pxq `∆zi grad gi pxq, ∆xy
iPE

iPA

“ xHessx Lpwqr∆xs, ∆xy `ÿ

∆yi xgrad hi pxq, ∆xy `iPEÿ∆zi xgrad gi pxq, ∆xy

iPA“ xHessx Lpwqr∆xs, ∆xy .
Thus, from second-order suﬃciency (A4), ∆x must be zero elements. And thenř ∆si “ 0 for all i P
IzA.
Next,
substituting
∆x
“
0
into
the
ﬁrst
equationin(11)yields0“
iPE ∆yi grad hi pxq `řiPA ∆zi grad gi pxq. The LICQ (A2) implies that the coeﬃcients ∆yi for i P E and ∆zi for i P A
must be zero. This completes the proof.

3.3

Prototype Algorithm of RIPM

Applying the Riemannian Newton method directly to the KKT vector ﬁeld F : N Ñ T N results inthe following Newton equation (see (3) without iteration count k) at each iteration:∇F pwqr∆ws ` F pwq “ 0.
5(12)

As with the usual interior point method in the Euclidean setting, once the iterates reach the boundary of the feasible region, they are forced to stick to it [33, P6]. For the iterates to maintain a
suﬃcient distance from the boundary, we introduce a perturbed complementary equation with somebarrier parameter µ ą 0 and deﬁne the perturbed KKT vector ﬁeld :Fµ pwq :“ F pwq ´ µê, and ê :“ êpwq :“ p0x , 0, 0, eq .

(13)

Notice that the perturbation term ê, indeed, is a special vector ﬁeld on N , not a constant, because
0x is essentially dependent on w and/or x. In fact, the covariant derivative of the perturbed KKT
vector ﬁeld is the same as that of the original. From the linearity of the connection ∇, for any
w P N and any µ ą 0, we have∇Fµ pwq “ ∇F pwq ´ µ∇êpwq “ ∇F pwq,

(14)

where the last equity comes from ∇êpwqr∆ws “ p0x , 0, 0, 0q for all ∆w P Tw N . Applying the
Riemannian Newton method to perturbed KKT vector ﬁeld Fµ pwq yields the perturbed Newton
equation: ∇Fµ pwqr∆ws ` Fµ pwq “ 0. From (13) and (14), this equation is equivalent to∇F pwqr∆ws ` F pwq “ µê,(15)which reduces to the ordinary Newton equation (12) as µ Ñ 0. At this point, we can describe aprototype of the Riemannian Interior Point Method (RIPM) in Algorithm 2.
Algorithm 2: Prototype Algorithm of RIPM for (RCOP)
Input: An initial point w0 “ px0 , y0 , z0 , s0 q P N with pz0 , s0 q ą 0 and a retraction R on M.
γ̂ P p0, 1q, µ0 ą 0.
Output: Sequence twk u Ă N such that twk u Ñ w˚ and w˚ satisﬁes the KKT conditions(4).
Set k Ñ 0;
while Stopping criterion not satisﬁed do
1. Obtain ∆wk “ p∆xk , ∆yk , ∆zk , ∆sk q P Twk N by solving the perturbed Newton
equation on Twk N :∇F pwk qr∆wk s “ ´F pwk q ` µk ê;(16)
2. Choose γk P rγ̂, 1s and compute the step size αk deﬁned by
"
"
*
"
**
psk qipzk qimin 1, γk min ´| p∆sk qi ă 0 , γk min ´| p∆zk qi ă 0;
iip∆sk qip∆zk qi(17)

3. Compute the next point as wk`1 :“ R̄wk pαk ∆wk q, see (6);
4. Choose 0 ă µk`1 ă µk ;5. k Ñ k ` 1;
end
In the step 2 of Algorithm 2, we just wish to compute a step size 0 ă αk ď 1 to ensure new
point wk`1 with pzk`1 , sk`1 q ą 0. There are many schemes to achieve this purpose. The scheme in(17) is simple but suﬃcient to guarantee the local convergence, as will be proved in Section 5.

3.4

Solving Perturbed Newton Equation Eﬃciently

The challenge of Algorithm 2 is how to solve the Newton equation (16) in an eﬃcient manner.
In this subsection, we will do this in two steps: the ﬁrst step will be to turn the original full
Newton equation, which is asymmetric and consists of four variables, into a condensed form, which
is symmetric and consists of only two variables. In the second step, an iterative method, namely,
the Krylov subspace method, is used to solve the operator equations directly, avoiding the expensive
computational eﬀort of converting them into the usual matrix equations.
3.4.1

Condensed Form of Perturbed Newton Equation

Let us consider Algorithm 2 and omit the iteration count k. Given the current point w P N withpz, sq ą 0, for the KKT vector ﬁeld F pwq in (5), we denote its components by Fx , Fy , Fz , Fs in topto-bottom order, namely, Fx :“ gradx Lpwq, Fy :“ hpxq, Fz :“ gpxq ` s, Fs :“ ZSe. By using these6

symbols and the formulation of ∇F pwq in (8), the full (perturbed) Newton equation (16) deﬁned
on Tw N – Tx M ˆ Rl ˆ Rm ˆ Rm is expanded as:¨˛ ¨˛Hessx Lpwqr∆xs ` Hx r∆ys ` Gx r∆zs
´Fx
˚ Hx˚ r∆xs
‹ ˚ ´Fy
‹
˚ ˚
‹“˚
‹.
(18)
˝ Gx r∆xs ` ∆s‚ ˝ ´Fz‚Z∆s ` S∆z´Fs ` µe
Not only does this equation contain four variables, but there is no symmetry on the left side ofthe equation, so it would be unwise to solve it just like that. By pz, sq ą 0 and the fourth line of(18), we deduce ∆s “ Z ´1 pµe ´ Fs ´ S∆zq . Substituting this ∆s into the third line of (18), we get
Gx˚ r∆xs´Z ´1 S∆z “ ´Z ´1 µe´gpxq and thus ∆z “ S ´1 rZ pGx˚ r∆xs ` Fz q ` µe ´ Fs s. Substituting
this ∆z further into the ﬁrst line of (18) and combining it with the second line of (18) yields thefollowing condensed Newton equation, which is deﬁned on Tx M ˆ Rl :ˆ˙ ˆ ˙
Aw r∆xs ` Hx r∆ysc
T p∆x, ∆yq :““.
(19)
Hx˚ r∆xsq
where

Aw :“ Hessx Lpwq ` Gx S ´1 ZGx˚ ,c :“ ´Fx ´ Gx S ´1 pZFz ` µe ´ Fs q , q :“ ´Fy .

(20)

Here, c and q are constant vectors. If we deﬁned Θ :“ Gx S ´1 ZGx˚ , then Aw “ Hessx Lpwq ` Θ. Notethat both Θ and Hessx Lpwq are operators from and to Tx M.
From the discussion above, for any w P N with pz, sq ą 0, the operator ∇F pwq in (8) is
nonsingular if and only if the newly deﬁned operator T in (19) is nonsingular. Eventually, it issuﬃcient for us to solve the equation (19) containing only two variables ∆x and ∆y. In fact, when
we consider the case of only inequality constraints in (RCOP), then ∆y vanishes, and only a linear
equation
Aw r∆xs “ c(21)
on Tx M needs to be solved. More importantly, the operator T in the left side of (19) is symmetric,or say self-adjoint (although often indeﬁnite). It is trivial to check that operators Θ and Aw are
self-adjoint on Tx M; and thus T is self-adjoint on the product vector space Tx M ˆ Rl equippedwith inner product xpξx , ξy q, pηx , ηy qy :“ xξx , ηx yx ` ξyT ηy . We can also see that (19) is a saddle point
problem deﬁned on Hilbert spaces form its special structure.
3.4.2

Krylov Subspace Methods on Tangent Space

Next, how to solve (19) eﬃciently becomes critical. For simplicity, we consider the case of only
inequality constraints in (RCOP), then we will solve operator equation (21) with a self-adjoint operator Aw : Tx M Ñ Tx M. Let d :“ dim Tx M. Unfortunately, in most cases of practical applications,
the Riemannian situation leaves us with no explicit matrix form available for Aw . This means that
we can only access A (subscript w omitted) by inputting a vector v to return Av.
A general approach is ﬁrst to ﬁnd the matrix representation Â for A under some basis of Tx M.
In detail, the full process of this approach is as follows:(Step 1) Obtain d random independent vectors on Tx M.
d(Step 2) Obtain an orthonormal basis tui ui“1 of Tx M by the modiﬁed Gram-Schmidt algorithm.
(Step 3) Compute rÂsij :“ xAuj , ui yx for 1 ď i ď j ď d due to symmetry, then we obtain thematrix representation Â P Rdˆd .
(Step 4) Compute rĉsi :“ xc, ui yx for 1 ď i ď d, then we obtain the vector representation ĉ P Rd .
(Step 5) Using arbitrary linear solver to get solution ∆x̂ P Rd from matrix equation Â∆x̂ “ ĉ.
řd(Step 6) Recovery the tangent vector ∆x P Tx M by ∆x “ i“1 p∆x̂qi ui .

7

In Algorithm 2, at each iteration, x is updated, and thus, the tangent space Tx M changes. Thus,
the above six steps need to be done all over again. Obviously, this approach is so expensive that it
is not feasible in practice.
An ideal approach is to use an iterative method, such as the Krylov subspace method (e.g.,
conjugate gradients method [5, Chapter 6.3]), on Tx M directly. Such a method does not explicitly
require a matrix representation Â for A. In general, it only needs to call an abstract linear operatorv ÞÑ Av. Since A in (20) is self-adjoint but indeﬁnite, for solving operator equation (21), we will
use the Conjugate Residual (CR) method (see [29, ALGORITHM 6.20]) as stated in Algorithm 3.
Algorithm 3: Conjugate Residual (CR) Method on Tangent Spaces for (21)Input: Symmetric invertible linear operator A : Tx M Ñ Tx M, nonzero c P Tx M and an
initial point v0 P Tx M.
Output: Sequence tvn u Ă Tx M such that tvn u Ñ v ˚ and Av ˚ “ c.
Set n Ñ 0, r0 :“ c ´ Av0 , p0 :“ r0 and compute Ar0 , Ap0 ;while stopping criterion not satisﬁed do
1. Update number αn :“ xrn , Arn yx {xApn , Apn yx ;// Step length
2. Set vn`1 :“ vn ` αn pn ;// Iterate point
3. Update rn`1 :“ rn ´ αn Apn ;// Residual
4. Compute Arn`1 ;// This is the only call to A in while loop
5. Update number βn :“ xrn`1 , Arn`1 yx {xrn , Arn yx ;// Conjugate direction
6. Set pn`1 :“ rn`1 ` βn pn ;7. Compute Apn`1 :“ Arn`1 ` βn Apn ;// No need to call A here
8. n Ñ n ` 1;
end
A signiﬁcant feature is that the iterate points vk , conjugate directions pk , and residual vectors
rk :“ Avk ´ c are all contained in Tx M. Usually, the initial point v0 is the zero element of Tx M;
the iteration terminates when the relative residual }rk } {}c} ď ϵ for some threshold ϵ ą 0, or somemaximum number of iterates is reached.
The discussion of the above two approaches can be naturally extended to the case containing
equality constraints in (RCOP), where we consider T on Tx M ˆ Rl instead of A on Tx M.

4Preliminaries and Auxiliary Results

This section introduces the useful results that are indispensable to our subsequent discussion.
Remark 4.1. For a retraction R on M and x P M, by DRx p0x q “ idTx M and the inverse function
theorem, there exists a neighborhood V of 0x in Tx M such that Rx is a diﬀeomorphism on V ; thus,
R´1
x pyq is well deﬁned for all y P M suﬃciently close to x. In this case, Rx pV q Ă M is called a
retractive neighborhood of x. Furthermore, the existence of a totally retractive neighborhood [38,Theorem 2] shows that for any x̄ P M there is a neighborhood W of x̄ such that R´1
x pyq is welldeﬁned for all x, y P W . In what follows, we will suppose that an appropriate neighborhood has
been chosen by default for the well-deﬁnedness of R´1
x pyq.

4.1

Vector Transport and Parallel Transport

Deﬁne the Whitney sum T M ‘ T M :“ tpη, ξq : η, ξ P Tx M, x P Mu . A smooth map T : T M ‘
T M Ñ T M : pη, ξq ÞÑ Tη pξq, is called a vector transport on M if there exists an associatedretraction R on M such that T satisﬁes the following properties for all x P M:
1. Associated retraction. Tη pξq P TRx pηq M for all η, ξ P Tx M.
2. Consistency. T0x pξq “ ξ for all ξ P Tx M.
3. Linearity. Tη paξ ` bζq “ aTη pξq ` bTη pζq for all a, b P R and η, ξ, ζ P Tx M.
Thus, ﬁxing any η P Tx M, the map Tη : Tx M Ñ TRx pηq M : ξ ÞÑ Tη pξq, is a linear operator.
Additionally, T is isometric if xTη pξq, Tη pζqy “ xξ, ζy holds, for all x P M and all ξ, ζ, η P Tx M. In
other words, for any η P Tx M, the adjoint and the inverse of Tη coincide, i.e., T˚η “ T´1
η . There
are two important classes of vector transport as follows. Let R be a retraction on M.
8

(1) The diﬀerentiated retraction deﬁned byTη pξq :“ DRx pηqrξs, η, ξ P Tx M, x P M,

(22)

is a valid vector transport [1, equation (8.6)].
(2) Given a smooth curve γ : r0, 1s Ñ M and t0 , t1 P r0, 1s, the parallel transport from the tangent
space at γpt0 q to the tangent space at γpt1 q along γ, is a linear operator Pγt1 Ñt0 : Tγpt0 q M Ñ
Tγpt1 q M deﬁned by Ptγ1 Ñt0 pξq “ Z pt1 q, where Z is the unique parallel vector ﬁeld such that
Z pt0 q “ ξ. Then, for any x P M, η P Tx M, thenTη pξq :“ P1Ñ0
pξq
γ

(23)

is a valid vector transport [1, equation (8.2)], where Pγ denotes the parallel transport along
the curve t ÞÑ γptq :“ Rx ptηq. We often omit the superscript 1Ñ0 if it is clear from the context.
In particular, parallel transport is isometric.

4.2

Lipschitz Continuity with respect to Vector Transports

Multiple Riemannian versions of Lipschitz continuity have been deﬁned, e.g., [5, Section 10.4]. Here,we consider Lipschitz continuity with respect to a vector transport. In what follows, let M be a
Riemannian manifold endowed with a vector transport T and an associated retraction R. We ﬁrst
consider the Lipschitz continuous gradient of scale ﬁeld f .
Deﬁnition 4.2 ([17, Deﬁnition 5.2.1]). A function f : M Ñ R is Lipschitz continuously diﬀerentiablewith respect to T in U Ă M if it is diﬀerentiable and there exists a constant κ ą 0 such that, for all
x, y P U , }grad f pyq ´ Tη rgrad f pxqs} ď κ}η}, where η “ R´1x y.
Going one degree higher, let us now discuss the Lipschitz continuity of Hessian operators.
Throughout this paper, for a linear operator A : E Ñ E 1 between two ﬁnite-dimensional normed
vector spaces E and E 1 , the (operator) norm of A is deﬁned by }A} :“ sup t}Av}E 1 : v P E, }v}E “ 1 ,or, }v}E ď 1u . The inverse of Tη is needed in the following deﬁnitions, so we can assume that vectortransport T is isometric, e.g., parallel transport in (23). In fact, there are many ways to constructisometric vector transports; see [18, Section 2.3].
Deﬁnition 4.3 ([18, Assumption 3]). A function f : M Ñ R is twice Lipschitz continuously diﬀerentiable with respect to T in U Ă M if it is twice diﬀerentiable and there exists a constant κ ą 0´1such that, for all x, y P U , } Hess f pyq ´ Tη Hess f pxqT´1η } ď κdpx, yq, where η “ Rx y.
Lemma 4.4 ([18, Lemma 4]). If f : M Ñ R is C 3 , then for any x̄ P M and any isometricvector transport T, there exists a neighborhood U of x̄ such that f is twice Lipschitz continuously
diﬀerentiable with respect to T in U .
If the operator, Hess f pxq, above is replaced by a general covariant derivative ∇F pxq, we can get
the next results in a similar way. Lemma 4.6 can be proven in the same way as Lemma 4.4.
Deﬁnition 4.5. Given a vector ﬁeld F on M. The map x ÞÑ ∇F pxq is Lipschitz continuous withrespect to T in U Ă M if there exists a constant κ ą 0 such that, for all x, y P U , it holds that´1}∇F pyq ´ Tη ∇F pxqT´1η } ď κdpx, yq, where η “ Rx y.
Lemma 4.6. If F is a C 2 vector ﬁeld, then for any x̄ P M and any isometric vector transport
T, there exists a neighborhood U of x̄ such that the map x ÞÑ ∇F pxq is Lipschitz continuous withrespect to T in U .

4.3

Auxiliary Lemmas

Notice that in the previous subsection on the deﬁnitions of Lipschitz continuity, we used }η} with
η “ R´1x y or dpx, yq to denote the upper bound on the right-hand side. The next lemma shows thatthe two are not essentially diﬀerent. When M “ Rn , both reduce to }x ´ y}.
Lemma 4.7 ([18, Lemma 2]). Let M be a Riemannian manifold with a retraction R and let x̄ P M.
Then,

9

(i) there exist a0 , a1 , δa0 ,a1 ą 0 such that for all x in a suﬃciently small neighborhood of x̄ andall ξ, η P Tx M with }ξ}, }η} ď δa0 ,a1 , one has a0 }ξ ´ η} ď dpRx pηq, Rx pξqq ď a1 }ξ ´ η}. In
particular, a0 }ξ} ď dpx, Rx pξqq ď a1 }ξ} when η “ 0;(ii) there exist a0 , a1 ą 0 such that for all x in a suﬃciently small neighborhood of x̄, one hasa0 }ξ} ď dpx, x̄q ď a1 }ξ} where ξ “ R´1x̄ pxq.
The next lemma is the fundamental theorem of calculus in the Riemannian case.
Lemma 4.8 ([18, Lemma 8]). Let F be a C 1 vector ﬁeld and x̄ P M. Then there exist a neighborhood
U of x̄ and a constant c1 ě 0 such that for all x, y P U ,
›
›
ż1
› 0Ñ1
›
0Ñt
tÑ0
›Pγ rF pyqs ´ F pxq ´› ď c1 }η}2 ,P∇F
pγptqqP
rηsdt
γ
γ
›
›
0

where η “ R´1x pyq and Pγ is the parallel transport along the curve γptq :“ Rx ptηq. Moreover, if
R “ Exp, then indeed c1 “ 0 above (see [13, equation (2.4)]).
The next lemma is a Riemannian extension of some important estimates, usually used to analyze
Newton methods [8, Lemma 4.1.12].
Lemma 4.9. Let F be a C 2 vector ﬁeld and x̄ P M. Then there exist a neighborhood U of x̄ and aconstant c2 ą 0 such that for all x P U ,› 0Ñ1››Pγ rF pxqs ´ F px̄q ´ ∇F px̄qrηs› ď c2 d2 px̄, xq,
where η “ R´1x̄ x and Pγ is the parallel transport along the curve γptq :“ Rx̄ ptηq.
Proof. Let LHS :“ }P0Ñ1
rF pxqs ´ F px̄q ´ ∇F px̄qrηs}. It follows that
γ
›
›
ż1
› 0Ñ1
›
0Ñt
tÑ0
›
LHS ď ›Pγ rF pxqs ´ F px̄q ´
Pγ ∇F pγptqqPγ rηsdt››
0
›
›ż 1
›
›
0Ñt
tÑ0
›
` › Pγ ∇F pγptqqPγ rηsdt ´ ∇F px̄qrηs››0›ż 1 `›˘››tÑ0ď c1 }η}2 ` ›P0Ñt∇FpγptqqP
´
∇F
px̄q
rηsdt
›. (by Lemma 4.8)
γ
γ
0

›
›ş `
˘
›
› 1
tÑ0
´
∇F
px̄q
rηsdt
∇F
pγptqqP
Let θ :“ › 0 P0Ñt›. Note that
γ
γ
ż1
θď

›
› 0Ñt
›Pγ ∇F pγptqqPtÑ0´ ∇F px̄q› }η} dt
γ0ż1c0 dpx̄, Rx̄ ptηqq }η} dt (by Lemma 4.6)

ď0
ż1

c0 a1 t}η} }η} dt “ď01
c0 a1 }η}2 . (by (i) of Lemma 4.7)
2

Combining the above results yields
1
1
LHS ď pc1 ` c0 a1 q}η}2 ď pc1 ` c0 a1 q{a20 d2 px̄, xq,
22
where the last inequality comes from (ii) of Lemma 4.7. Letting c2 :“ pc1 ` 12 c0 a1 q{a20 completes
the proof.
We end this section with the following useful lemmas.
Lemma 4.10 ([12, Lemma 3.2]). Given a vector ﬁeld F on M. If the map p ÞÑ ∇F ppq is continuousat p˚ and ∇F pp˚ q is nonsingular, then there exist›a neighborhood
U of p˚ and a constant Ξ ą 0›´1such that, for all p P U , ∇F ppq is nonsingular and ›∇F ppq › ď Ξ.
Lemma 4.11 ([14, Lemma 14.5]). Let F be a C 2 vector ﬁeld on M and p˚ P M. If F pp˚ q “ 0 and∇F pp˚ q is nonsingular, then there exist a neighborhood U of p˚ and constants c3 , c4 ą 0 such that,for all p P U , c3 dpp, p˚ q ď }F ppq} ď c4 dpp, p˚ q.
Lemma 4.12 ([7, Lemma 3.5]). Let u P Tp M such that Exp› p puq exists and› v P Tp M – Tp pTp Mq.
Then xD Expp puqrus, D Expp puqrvsy “ xu, vy. In particular, ›D Expp pλuqrus› “ }u} holds all λ ě 0.
10

5Local Convergence

Here, for any two nonnegative sequences tuk u and tvk u, we write uk “ Opvk q if there is a constantM ą 0 such that uk ď M vk for all suﬃciently large k; and we write uk “ opvk q if vk ą 0 and thesequence of ratios tuk {vk u approaches zero. In this section, we will establish local convergence ofour prototype Algorithm 2 of RIPM.

5.1

Perturbed Damped Riemannian Newton Method

We will rely on an application of the so-called perturbed damped Riemannian Newton method forsolving the singularity problem (2), which can be stated as Algorithm 4.
Algorithm 4: Perturbed Damped Riemannian Newton Method for (2)Input: A vector ﬁeld F P XpMq, an initial point p0 P M and a retraction R on M. Setµ0 ą 0.
Output: Sequence tpk u Ă M such that tpk u Ñ p˚ and F pp˚ q “ 0p˚ .
Set k Ñ 0;
while Stopping criterion not satisﬁed do
1. Obtain ξk P Tpk M by solving the perturbed Newton equation:∇F ppk qrξk s “ ´F ppk q ` µk ê;(24)

2. Choose a (damped) step size 0 ă αk ď 1;
3. Compute the next point as pk`1 :“ Rpk pαk ξk q;
4. Choose 0 ă µk`1 ă µk ;5. k Ñ k ` 1;
end
In contrast to the standard Riemannian Newton method described in Algorithm 1, the term
“perturbed” means that we solve a Newton equation with a perturbed term µk ê, while “damped”
means using αk instead of unit steps. It is well known that Algorithm 1 are locally superlinearly[12] and quadratically [13] convergent under the following Riemannian Newton assumptions:(B1) There exists p˚ P M such that F pp˚ q “ 0p˚ ;(B2) The covariant derivative ∇F pp˚ q is nonsingular;(B3) The vector ﬁeld F is C 2 .
As the following Proposition 5.1 shows, Algorithm 4 also has the same convergence properties
as Algorithm 1 if we control µk and αk according to the two schemes that Proposition 5.1 gives.
We can see that either scheme will have µk Ñ 0 and αk Ñ 1, which makes Algorithm 4 eventually
reduce to Algorithm 1 when k is suﬃciently large.
Proposition 5.1 (Local convergence of Algorithm 4). Consider the perturbed damped Riemannian
Newton method described in Algorithm 4 for the singularity problem (2). Let (B1)-(B3) hold.
Choose parameters µk , αk as follows; then there exists a constant δ ą 0 such that for all p0 P M
with dpp0 , p˚ q ă δ, the sequence tpk u is well deﬁned. Furthermore,(1) if we choose µk “ op}F ppk q}q and αk Ñ 1, then pk Ñ p˚ superlinearly;(2) if we choose µk “ Op}F ppk q}2 q and 1 ´ αk “ Op}F ppk q}q, then pk Ñ p˚ quadratically.
˚
Proof. By (B2)-(B3), Lemma 4.6› and Lemma› 4.10, we can let pk be suﬃciently close to p such
´1 ››
ď Ξ for some constant Ξ. Then, the next iteratethat ∇F ppk q is nonsingular, and ∇F ppk q
point,
pk`1 :“ Rpk rαk ∇F ppk q´1 p´F ppk q ` µk êqs,
˚
is well deﬁned in Algorithm 4, and it follows from p˚ “ Rpk pηq with η :“ R´1
pk p and (i) of Lemma
4.7 thatdppk`1 , p˚ q ď a1 }η ´ αk ∇F ppk q´1 p´F ppk q ` µk êq}
“ a1 }η ` αk ∇F ppk q´1 pF ppk q ´ µk êq}.
11

(25)

Let rk :“ η ` αk ∇F ppk q´1 pF ppk q ´ µk êq. Algebraic manipulations show thatrk “ p1 ´ αk qη ` αk ∇F ppk q´1 r∇F ppk qη ` F ppk q ´ P0Ñ1
F pp˚ q ´ µk ês,γwhere Pγ is the parallel transport along the curve γptq “ Rpk ptηq and F pp˚ q “ 0. Thus, using}η} ď dppk , p˚ q{a0 from (ii) of Lemma 4.7 and Lemma 4.9, we have}rk } ď p1 ´ αk q}η} ` αk }∇F ppk q´1 }}P0Ñ1
F pp˚ q ´ F ppk q ´ ∇F ppk qη}γ` αk }∇F ppk q´1 }}ê}µk
ď p1 ´ αk qdppk , p˚ q{a0 ` αk }∇F ppk q´1 }c2 d2 ppk , p˚ q
` αk }∇F ppk q´1 }}ê}µk
ď p1 ´ αk qdppk , p˚ q{a0 ` Ξc2 d2 ppk , p˚ q ` Ξ}ê}µk .
Combining the above with (25), we conclude thatdppk`1 , p˚ q ď κ1 p1 ´ αk qdppk , p˚ q ` κ2 d2 ppk , p˚ q ` κ3 µk(26)

for some positive constants κ1 , κ2 , κ3 . On the other hand, by Lemma 4.11, we have
}F ppk q} “ Opdppk , p˚ qq.

(27)

In what follows, we prove assertions (1) and (2).
(1) Suppose that αk Ñ 1 and µk “ op}F ppk q}q, which together with (27) imply µk “ opdppk , p˚ qq.
By (26), we havedppk`1 , p˚ qµk
ď κ1 p1 ´ αk q ` κ2 dppk , p˚ q ` κ3.
(28)
dppk , p˚ q
dppk , p˚ q
We can take δ suﬃciently small and k suﬃciently large, if necessary, to conclude that dppk`1 , p˚ q ă1˚˚˚2 dppk , p q ă δ. Thus, pk`1 P Bδ pp q, the open ball of radius δ centered at p on M. By induction,
it is easy to show that the sequence tpk u is well deﬁned and converges to p˚ . Taking the limit ofboth sides of (28) proves superlinear convergence.
(2) Again, we start from (26) and rewrite it as:
dppk`1 , p˚ q “ p1 ´ αk qOpdppk , p˚ qq ` Opd2 ppk , p˚ qq ` Opµk q.

(29)

2

Suppose that 1 ´ αk “ Op}F ppk q}q and µk “ Op}F ppk q} q. Using (27), the above reduces todppk`1 , p˚ q “ Opd2 ppk , p˚ qq. This implies that there exists a constant ν such that dppk`1 , p˚ q ď
νd2 ppk , p˚ q, and hence, dppk`1 , p˚ q ď νd2 ppk , p˚ q ď νδ 2 ă δ, if δ is suﬃciently small. Again, by
induction, tpk u converges to p˚ quadratically.

5.2

Local Convergence of Algorithm 2

Next, lemma shows the relationship between the parameter γk and step size αk in Algorithm 2.
Lemma 5.2. Consider the Algorithm 2 for solving the problem (RCOP). Let (A1) and (A3) holdat some w˚ “ px˚ , y ˚ , z ˚ , s˚ q and αk be as in (17). Deﬁne a constant,!
)Π :“ 2 max max t1{ps˚ qi | ps˚ qi ą 0u , max t1{pz ˚ qi | pz ˚ qi ą 0u .
ii

For γk P p0, 1q, if Π }∆wk } ď γk , then 0 ď 1 ´ αk ď p1 ´ γk q ` Π }∆wk } .
Proof. Notice that the fourth line of (18) yields
Sk´1 ∆sk ` Zk´1 ∆zk “ µk pSk Zk q´1 e ´ e,which is exactly the same as in the usual interior point method in the Euclidean setting. Thus,
the proof entails directly applying [35, Lemma 3 and 4] for the Euclidean case to the Riemanniancase.
Now, let us establish the local convergence of our Algorithm 2 in a way that replicates Proposition
5.1 except for taking account of parameter γk .
12

Theorem 5.3 (Local convergence of prototype Algorithm 2). Consider the Algorithm 2 for solving
the problem (RCOP). Let (A1)-(A4) hold at some w˚ . Choose parameters µk , γk as follows; thenthere exists a constant δ ą 0 such that, for all w0 P N with dpw0 , w˚ q ă δ, the sequence twk u is welldeﬁned. Furthermore,(1) if we choose µk “ op}F pwk q}q and γk Ñ 1, then wk Ñ w˚ superlinearly;
2(2) if we choose µk “ Op}F pwk q} q and 1 ´ γk “ Op}F pwk q}q, then wk Ñ w˚ quadratically.
Proof. We only prove (2) because (1) can be proven in the same way. Let wk be such thatdpwk , w˚ q ă δ for suﬃciently small δ. From Proposition 3.1, (A1)-(A4) shows that the KKT vector ﬁeld F satisﬁes (B1)-(B3); thus, the discussion in the proof of Proposition 5.1 applies to KKT
vector ﬁeld F as well, simply by replacing the symbol p with w. Since we choose µk “ Op}F pwk q}2 q,and }F pwk q} “ Opdpwk , w˚ qq by (27), we obtain µk “ Opd2 pwk , w˚ qq. Thus,››}∆wk } “ ›∇F pwk q´1 p´F pwk q ` µk êq› (by (16))ď Ξp}F pwk q} ` µk }ê}q (by Lemma 4.10)
ď Op}F pwk q}q ` Opµk q“ Opdpwk , w˚ qq ` Opd2 pwk , w˚ qq “ Opdpwk , w˚ qq.
Since δ is suﬃciently small, from the above inequalities, the condition of Lemma 5.2 is satisﬁed,
i.e., }∆wk } ď γ̂{Π ď γk {Π for a constant γ̂ P p0, 1q. Hence, by 1 ´ γk “ Op}F pwk q}q, one has
1 ´ αk ď p1 ´ γk q ` Π }∆wk } “ p1 ´ γk q ` Opdpwk , w˚ qq “ Opdpwk , w˚ qq. Finally, from (29), we have
dpwk`1 , w˚ q “ p1 ´ αk qOpdpwk , w˚ qq ` Opd2 pwk , w˚ qq ` Opµk q “ Opd2 pwk , w˚ qq. This completes theproof.
Algorithm 2 guaranteed the local convergence, but we are more interested in its globally convergent version (described in the next section). We still provide a simple example of Algorithm 2 online22, where from (2) of Theorem 5.3, in practice we set parameters µk`1 “ mintµk {1.5, 0.5 }F pwk q} u;and γk`1 “ maxtγ̂, 1 ´ }F pwk q}u with γ̂ “ 0.5.

6

Global Algorithm

The globally convergent version of our RIPM uses the classical line search described in [11]. The
following considerations and deﬁnitions are needed in order to describe it compactly. For simplicity,
we often omit the subscript of iteration count k.
Given the current point w “ px, y, z, sq and ∆w “ p∆x, ∆y, ∆z, ∆sq, the next iterate is obtained
along a curve on product manifold N , i.e., α ÞÑ wpαq :“ R̄w pα∆wq with some step size α ą 0, see (6)for R̄w . By introducing wpαq “ pxpαq, ypαq, zpαq, spαqq, we have xpαq “ Rx pα∆xq, ypαq “ y ` α∆y,zpαq “ z ` α∆z, and spαq “ s ` α∆s. For a given starting point w0 P N with pz0 , s0 q ą 0, let us
set two constants
τ1 :“ minpZ0 S0 eq{pz0T s0 {mq, τ2 :“ z0T s0 { }F pw0 q} .
As well, deﬁne two functions f I pαq :“ minpZpαqSpαqeq ´ γτ1 zpαqT spαq{m, f II pαq :“ zpαqT spαq ´γτ2 }F pwpαqq}, where γ P p0, 1q is a constant. For i “ I, II, deﬁne␣(αi :“ max α : f i ptq ě 0, for all t P p0, αs ,(30)
αPp0,1si.e., αi are either one or the smallest positive root for the functions f i pαq in p0, 1s.
Moreover, we deﬁne the merit function φ : N Ñ R by φpwq :“ }F pwq}2 ; accordingly, we have
grad φpwq “ 2∇F pwq˚ rF pwqs, where symbol ˚ means its adjoint operator. Let } ¨ }1 , } ¨ }2 be l1 , l2222vector norms. Note that }F pwq}2w “ }gradx Lpwq}x `}hpxq}22 `}gpxq`s}? 2 `}ZSe}2 by (5). Moreover,mT
for any nonnegative z, s P R , one has }ZSe}2 ď z s “ }ZSe}1 ď m }ZSe}2 . Hence,?
?
}ZSe}2 { m ď z T s{ m ď }ZSe}2 ď }F pwq}.
(31)
Now, using the above deﬁnitions, the globally convergent RIPM can be stated as Algorithm 5.
2 See .../LocalRIPM/PrototypeRIPM.m in https://doi.org/10.5281/zenodo.10612799

13

Algorithm 5: Global Convergent Algorithm of RIPM for (RCOP)
Input: An initial point w0 “ px0 , y0 , z0 , s0 q P N with pz0 , s0 q ą 0 and a retraction R on M.
θ P p0, 1q, β P p0, 0.5s, γ´1 P p0.5, 1q.
Output: Sequence twk u Ă N such that twk u Ñ w˚ and w˚ satisﬁes the KKT conditions(4).
Set k Ñ 0;
while Stopping criterion not satisﬁed do?
1. Set σk P p0, 1q, ρk P rzkT sk {m, }F pwk q} { ms;2. Obtain ∆wk “ p∆xk , ∆yk , ∆zk , ∆sk q P Twk N by solving the following linear equation:∇F pwk qr∆wk s “ ´F pwk q ` σk ρk ê;(32)

3. Step size selection:(3a) Centrality condition: Set γk P p0.5, γk´1 q and ᾱk “ mintαkI , αkII u from (30);(3b) Suﬃcient decrease condition: Let αk :“ ᾱk ;while αk does not satisfy the condition φpR̄wk pαk ∆wk qq ´ φpwk q ď αk β xgrad φk , ∆wk ydo
αk :“ θαk ;end
4. Compute the next point as wk`1 :“ R̄wk pαk ∆wk q;
5. k Ñ k ` 1;
end

Compared to the prototype Algorithm 2, the global Algorithm 5 involves a more elaboratechoice of step size. Regarding the centrality condition (3a), it does not diﬀer in any way from the
Euclidean setting; thus, references [11, 3, 9] show that ᾱk is well-deﬁned thereby ensuring that zk ,sk are positive. In the following, we focus on the suﬃcient decrease (Armijo) condition (3b):φpR̄wk pαk ∆wk qq ´ φpwk q ď αk β xgrad φk , ∆wk y ,(33)

where grad φk ” grad φ pwk q for short. With a slight abuse of notation` φ, at the
˘ current point w
and direction ∆w, we deﬁne a real-to-real function α ÞÑ φpαq :“ φ R̄w pα∆wq ; then, it follows
from the deﬁnition of a retraction and the chain rule that“‰
φ1 p0q “ DφpR̄w p0qq DR̄w p0qr∆ws “ Dφpwqr∆ws “ xgrad φpwq, ∆wy.
Hence, φ1k p0q “ xgrad φk , ∆wk y at the k-th iterate wk ; and then the Armijo condition (33) reads
φk pαk q ´ φk p0q ď αk βφ1k p0q as usual. If φ1k p0q ă 0, the backtracking loop in (3b) of Algorithm 5
will eventually stop [26, Lemma 3.1]. The next lemma shows the condition under which the Newton
direction ∆wk generated by (32) ensures the descent of the merit function.
Lemma 6.1. If the direction ∆wk is the solution of equation (32), thenxgrad φpwk q, ∆wk y “ 2p´}F pwk q}2 ` σk ρk zkT sk q.
In this case, ∆wk is a descent direction for φpwq at wk if and only if ρk ă }F pwk q}2 {σk zkT sk .
Proof. The iteration count k is omitted. Let ∆w be given by (32). Then, we have xgrad φpwq, ∆wy “x2∇F pwq˚ rF pwqs, ∆wy “ 2xF pwq, ∇F pwqr∆wsy “ 2xF pwq, ´F pwq ` σρêy “ 2p´xF pwq, F pwqy `σρxF pwq, êyq. Then, by deﬁnition of ê in (13), xF pwq, êy “ xZSe, ey “ z T s completes the proof.
The next proposition shows that Algorithm 5 can generate the monotonically nonincreasing
sequence tφk u. Note that φk`1 ” φk pαk q and φk ” φk p0q ” φ pwk q.
Proposition 6.2. If }F pwk q} ‰ 0, then the direction ∆wk generated by Algorithm 5 is a descentdirection for φpwq at wk . Moreover, if the Armijo condition (33) is satisﬁed, then
φk pαk q ď r1 ´ 2αk βp1 ´ σk qs φk p0q.
Thus, the sequence tφk u is monotonically nonincreasing.
14

?
Proof. The iteration count k is omitted. Suppose that ρ ď }F pwq} { m and ∆w is given by (32),
we have
φ1 p0q “ xgrad φpwq, ∆wy “ 2p´φpwq ` σρz T sq (by Lemma 6.1)?
ď 2p´φpwq ` σ }F pwq} z T s{ mq2ď 2p´φpwq ` σ }F pwq} q (by (31))“ ´2p1 ´ σqφpwq ă 0.

(34)

Thus, in Algorithm 5, ∆w is a descent direction
function φ at w. Alternatively,by? for the merit
?
2
T
T
T
Lemma 6.1,
itis
suﬃcient
to
show
that
}F
pwq}{
m
ă
}F
pwq}
{σz
s.
By
σz
s
ă
z
s
ď
m}F
pwq};
?
then, 1{ m ă }F pwq}{σz T s. Multiplying both sides by }F pwq} yields the result.
Moreover, if condition (33) is satisﬁed, then by (34), we have φpαq ď φp0q`αβxgrad φpwq, ∆wy ď
φp0q ` αβp´2p1 ´ σqφp0qq “ r1 ´ 2αβp1 ´ σqsφp0q. Note that in Algorithm 5, we set β P p0, 1{2s, σ P
p0, 1q, and α P p0, 1s, which imply that the sequence tφk u is monotonically nonincreasing.
Finally, we need to make the following assumptions. For ϵ ě 0, let␣(Ωpϵq :“ w P N | ϵ ď φpwq ď φ0 , minpZSeq{pz T s{mq ě τ1 {2, z T s{}F pwq} ě τ2 {2 .
l(C1) In the set Ωp0q, f, h, and g are smooth functions; tgrad hi pxqui“1 is linearly independent for
all x; and the map w ÞÑ ∇F pwq is Lipschitz continuous (with respect to parallel transport);(C2) The sequences txk u and tzk u are bounded [10, 3];(C3) In any compact subset of Ωp0q, ∇F pwq is nonsingular.
Given the above assumptions, we can now prove the following statement.
Theorem 6.3 (Global Convergence of RIPM). Let twk u be generated by Algorithm 5 with R “ Expand tσk u Ă p0, 1q be bounded away from zero and one. Let φ be Lipschitz continuous on Ωp0q. If(C1)-(C3) hold, then t}F pwk q}u converges to zero. Moreover, if w˚ is a limit point of sequence
twk u, then w˚ satisﬁes Riemannian KKT conditions (4).
The proof of the above theorem will be given in the next section. Note that although theexponential map is used in the proof, the numerical experiments indicate that global convergence
may hold for a general retraction R.

7

Proof of Global Convergence

In this section, our goal is to prove the global convergence Theorem 6.3. We will follow the proof
procedure in [11], which discussed Algorithm 5 when M “ Rn in (RCOP), and we call the algorithmin [11] the Euclidean Interior Point Method (EIPM). In what follows, we will omit similar content
because of space limitations and focus on the diﬃculties encountered when adapting the proof of
EIPM to RIPM. In particular, we will make these diﬃculties as tractable as in EIPM by proving aseries of propositions in Subsection 7.1.

7.1

Continuity of Some Special Scalar Fields

To show the boundedness of the sequences generated by Algorithm 5, we need the continuity ofsome special scalar ﬁelds on manifold M. The claims of this subsection are trivial if M “ Rn , butthey need to be treated carefully for general M.
If we assign a linear operator Ax : Tx M Ñ Tx M to each x P M, then the map x ÞÑ }Ax } :“sup t}Ax v}x | v P Tx M, }v}x “ 1, or, }v}x ď 1u is a scalar ﬁeld on M; but notice that the operator
norm } ¨ } depends on x. Let Sympdq denote the set of symmetric matrices of order d, and } ¨ }F ,} ¨ }2 denote the Frobenius norm and the spectral norm, respectively, applied to a given matrix.
Lemma 7.1. Let M be an n-dimensional Riemannian manifold. Let x P M and Ax be a linear
operator on Tx M. Choose an orthonormal basis of Tx M with respect to x¨, ¨yx , and let Âx P Rnˆndenote the matrix representation of Ax under the basis. Then, }Âx }2 , }Âx }F are invariant under achange of orthonormal basis; moreover, }Ax } “ }Âx }2 ď }Âx }F .
15

Proof. Suppose that there are two orthonormal bases tEi uni“1 , tEi1 uni“1 on Tx M. With respect tothem, let P P Rnˆn denote the change-of-basis matrix, i.e., rP skj :“ xEj1 , Ek yx , for 1 ď k, j ď n;then, P is orthogonal. Let Âx , Â1x P Rnˆn denote the matrix representations of Ax under thetwo bases, respectively. We have Â1x “ P ´1 Âx P. Then, }Â1x } “ }P ´1 Âx P } “ }Âx } holds for the
Frobenius norm or the spectral norm. Therefore, the values }Âx }2 and }Âx }F are invariant undera change of orthonormal basis. Now, consider an orthonormal řbasis tEi uni“1 on Tx M. For anynny P Tx M, its vector representation ŷ P R is deﬁned by y “ i“1 ŷi Ei . Accordingly, we haveřnyAx y “ Âx ŷ, i.e., Ax y “i“1 pÂx ŷqi Ei , and from the orthonormal property of the basis, we have2

}Ax y}x “ }nÿpÂx ŷqi Ei }2x “

i“1

n
ÿ

pÂx ŷq2i “ }Âx ŷ}22 .

i“1Thus, }Ax y}x “ }Âx ŷ}2 for any y P Tx M. Finally, we have}Ax } “supyPTx M,}y}x “1}Ax y}x “supŷPRn ,}ŷ}2 “1}Âx ŷ}2 “ }Âx }2 .

It is clear that }X}2 ď }X}F for any matrix X.
Yet, we have not clariﬁed the continuity about x ÞÑ }Ax }. The following proposition proves thecontinuity of an important case of Ax that appears in (RCOP).
ˆ pxq P Sympdq denote the matrix representation
Proposition 7.2. Consider f in (RCOP). Let Hessf
ˆ pxq} is a continuousof Hess f pxq under an arbitrary orthonormal basis of Tx M. Then, x ÞÑ }Hessf
scalar ﬁeld on M for } ¨ }F or } ¨ }2 . Moreover, x ÞÑ } Hess f pxq} is a continuous scalar ﬁeld on M.
ˆ pxq} is well deﬁned, so it suﬃces to prove its continuity.
Proof. Lemma 7.1 shows that x ÞÑ }Hessf
From Corollary 13.8 in [21], for each x̄ P M there is a smooth, orthonormal local frame tEi udi“1on a neighborhood U of x̄; namely, tE1 pxq, . . . , Ed pxqu forms an orthonormal basis on Tx M for allx P U . Choose such a local frame tEi udi“1 around x̄; then, the matrix representation of Hess f pxq isˆ pxqskj :“ xHess f pxqrEj pxqs, Ek pxqy “ xp∇E grad f qpxq, Ek pxqyx for 1 ď k, j ď d.
given by rHessfjxˆ pxq is a continuous
From the smoothness of the Riemannian metric (1), it follows that x ÞÑ Hessfˆfunction from U Ă M to Sympdq. Since matrix norms are continuous, }Hessf pxq} is continuous onˆ pxq}2 “ } Hess f pxq} for any
U Q x̄. This argument holds for any x̄ P M. From Lemma 7.1, }Hessf
x P M, which completes the proof.
lm

The above result can be applied verbatim to the Hessian of constraint functions thi ui“1 , tgi ui“1in (RCOP). The next proposition can be proved similarly, as in Lemma 7.1 and Proposition 7.2.
Proposition 7.3. Consider h, g in (RCOP) and the linear operators Hx , Gx deﬁned in (7). Then,x ÞÑ }Hx } and x ÞÑ }Gx } are continuous scalar ﬁelds on M.
Proposition 7.4. Given w “ px, y, z, sq P N , consider the operator ∇F pwq in (8). Let tEi udi“1 belman orthonormal basis of Tx M and tei uli“1 , te9 i um
i“1 be the standard bases of R , R , respectively. If
we choose an orthonormal basis of Tw N as follows:
d

l

m

m

tpEi , 0, 0, 0qui“1 Y tp0x , ei , 0, 0qui“1 Y tp0x , 0, e9 i , 0qui“1 Y tp0x , 0, 0, e9 i qui“1 ,then, the matrix representation of ∇F pwq is given by¨
Q B
˚ BT 0ˆ pwq “ ˚ T
∇F
˝ C
0
0
0

C
0
0
S

(35)

˛
0
0 ‹
‹,
I ‚
Z

i.e., a matrix of order pd ` l ` 2mq and, where
Q :“ Qpwq P Sympdq is given by rQskj :“ xHessx LpwqrEj s, Ek yx for 1 ď k, j ď d;
ˆ 1 pxq, ¨ ¨ ¨ , gradh
ˆ l pxqs P Rdˆl ;B :“ Bpxq “ rgradh
ˆ 1 pxq, ¨ ¨ ¨ , gradg
ˆ m pxqs P Rdˆm and the “hat” symbol above means the correC :“ Cpxq “ rgradg
sponding vector representation under the basis tEi udi“1 .
In this case, there is a continuous scalar ﬁeld T : N Ñ R such that for any w, }Qpwq}F ď T pwq.
Moreover, x ÞÑ }Bpxq}F and x ÞÑ }Cpxq}F are continuous scalar ﬁelds on M.
16ˆ pwq under the basis (35) is obtained through a trivial process, so we will omit
Proof. The matrix ∇Fˆ x Lpwq “ Hessf
ˆ pxq ` řl yi Hessh
ˆ i pxq `its description. From relation (9), we have Qpwq :“ Hessi“1řmˆi“1 zi Hessgi pxq, under the same basis. Thus,ˆ pxq}F `}Qpwq}F ď }Hessflÿˆ i pxq}F `|yi |}Hessh

i“1

m
ÿ

ˆ i pxq}F “: T pwq.
|zi |}Hessg

i“1ˆ pxq}F , }Hessh
ˆ i pxq}F , and }Hessg
ˆ i pxq}F are all continuous with respect
From Proposition 7.2, }Hessf
to x, thus, T is continuous. As for x ÞÑ }Bpxq}F , since the basis tEi udi“1 is orthonormal,}Bpxq}2F “lÿˆ i pxq}2 “}gradh
2

i“1

l
ÿ

2

}grad hi pxq}x ,i“1which implies continuity by (1). The claim for x ÞÑ }Cpxq}F can be proven similarly.

7.2

Global Convergence Theorem

Now, we are ready to prove the global convergence Theorem 6.3 by following the procedure in [11].
In what follows, we will omit similar content in [11] and focus on the diﬃculties encountered when
adapting the proof of EIPM to RIPM.
Proposition 7.5 (Boundedness of the sequences). Let twk u be a sequence generated by Algorithm
5 and suppose that (C1)-(C3) hold. If ϵ ą 0 and wk P Ωpϵq for all k, then(a) tzkT sk u, tpzk qi psk qi u , i “ 1, . . . , m, are all bounded above and below away from zero;(b) tzk u and tsk u are bounded above and component-wise bounded away from zero;(c) twk u is bounded;(d) t}∇F pwk q´1 }u is bounded;(e) t∆wk u is bounded.
Proof. The proofs in [11, Lemma 6.1] and/or [3, Theorem 2 (a)] can be applied verbatim to (a),(b) and (e).
(c) By (b), it suﬃces to prove that tyk u is bounded. The iteration count k is omitted in what
follows. By using the notation Hx and Gx as deﬁned in (7), we have Hx y “ gradx Lpwq´grad f pxq´
Gx z “: b. By (C1), Hx is an injection; then, there exists a unique solution y to Hx y “ b. Indeed,
we have´1y “ rpHx˚ Hx q Hx˚ s pgradx Lpwq ´ grad f pxq ´ Gx zq .
(36)
´1

Deﬁne Cx : Tx M Ñ Rl as Cx :“ pHx˚ Hx q Hx˚ . Under an arbitrary orthonormal basis of Tx M andstandard basis of Rl , if Ĥx is the matrix corresponding to Hx , then Cˆx “ pĤxT Ĥx q´1 ĤxT . It is easy
to show that }Cx } “ }Cˆx }2 for any x. We can see that for each x̄ P M there is a neighborhood U of x̄
such that x ÞÑ Ĥx is continuous over U . Then, by function composition, x ÞÑ Cˆx is also continuous
over U . This implies that x ÞÑ }Cx } “ }Cˆx }2 is continuous at each x̄, and hence, on M. Finally, by
Proposition 7.3, }Cx } , }grad f pxq} , and }Gx } are all continuous on M. Because`?txk u is bounded, by
˘(36) we have }yk } ď }Cxk } p}gradx Lpwk q} ` }grad f pxk q} ` }Gxk } }zk }q ď c1
φ0 ` c2 ` c3 }zk } ,for some positive constants c1 , c2 , c3 . Then, tyk u is bounded because tzk u is bounded.
(d) For each wk , choose an arbitrary orthonormal basis of Twk N . If the matrix representation
ˆ pwk q corresponds to ∇F pwk q, then r∇Fˆ pwk qs´1 corresponds to ∇F pwk q´1 . By Lemma 7.1, we
∇F
›
›
´1 ›
´1
ˆ
ˆ pwk qs´1 }F u is bounded.
›have ∇F pwk qď }r∇F pwk qs }F ; thus, it is suﬃcient to show that t}r∇F
For convenience, we will choose the basis of Twk N given in (35). Then, we have˛¨
Qk Bk Ck 0
˚ T
0
0
0 ‹
ˆ pwk q “ ˚ BkT
‹.
∇F
˝ Ck
0
0
I ‚
0
0 Sk Z k
By Proposition 7.4, there is a continuous scalar ﬁeld T : N Ñ R such that }Qpwq}F ď T pwq
for all w P N ; and }Bpxq}F , }Cpxq}F are continuous on M. It follows from the boundedness of17

txk u and twk u that for all k, }Qk }F ” }Qpwk q}F ď T pwk q ď c4 , }Bk }F ” }Bpxk q}F ď c5 , and}Ck }F ” }Cpxk q}F ď c6 , for some positive constants c4 , c5 , and c6 .
ˆ pwk q and theOn the other hand, whichever basis is used in the form of (35), the structure of ∇F
properties of its block submatrices remain unchanged, e.g., symmetry of Qk ; full rank of Bk ; identity
matrix I in the third row; all zero matrices; diagonal matrices Sk , Zk ; etc. This ensures that we canˆ pwk q. Up to this point,
obtain the desired result by performing an appropriate decomposition of ∇F
we have created all the conditions needed in the proof of the Euclidean version, namely, EIPM. Weˆ pwk qs´1 }F u is bounded by applying the proofs in [11, Lemma 6.2]
can make the claim that t}r∇F
and/or [3, Theorem 2 (c)] directly.
Lemma 7.6 (tᾱk u bounded away from zero). Let twk u be generated by Algorithm 5 with R “ Expand let (C1)-(C3) hold. If ϵ ą 0 and wk P Ωpϵq for all k, tσk u is bounded away from zero; then,tᾱk u is bounded away from zero.
Proof. Since ᾱk “ mintαkI , αkII u, it is suﬃcient to show that tαkI u and tαkII u are bounded away fromzero. For αkI , see [11, Lemma 6.3] and/or [9, Theorem 3.1]. The proofs in those references apply
verbatim to the Riemannian case. On the other hand, for αkII , we need to adapt the proofs in [11, 9],since Lipschitz continuity on manifolds is more complicated, see Subsection 4.2.
¯ w pα∆wq. Fix α∆w and let Pγ be theLet us suppress the subscript k. Recall that wpαq “ Exp¯parallel transport along the geodesic γptq “ Expw ptα∆wq. By Lemma 4.8 where c1 “ 0, we obtainPγ0Ñ1 rF pwpαqqs “ F pwq ` α∇F pwqr∆ws ´ α∇F pwqr∆ws `

ż1

P0Ñt
∇F pγptqqPtÑ0rα∆wsdtγγ

0

ż1“ F pwq ` α pσρê ´ F pwqq ` α` 0Ñt˘Pγ ∇F pγptqqPtÑ0´ ∇F pwq r∆wsdt
γ0

ż1“ p1 ´ αqF pwq ` ασρê ` α` 0Ñt˘Pγ ∇F pγptqqPtÑ0´ ∇F pwq r∆wsdt.
γ

0Taking the norm on both sides above gives››}F pwpαqq} “ ›P0Ñ1
rF pwpαqqs› (since parallel transport (23) is isometric)
γ
ż1
› 0Ñt
›
›Pγ ∇F pγptqqPtÑ0ď p1 ´ αq }F pwq} ` ασρ }ê} ` α´ ∇F pwq› }∆w} dt
γ
0

?

ż1

ď p1 ´ αq }F pwq} ` ασρ m ` α
κ2 }tα∆w} }∆w} dt0?
“ p1 ´ αq }F pwq} ` ασρ m ` α2 }∆w}2 κ{2.
The rest of the proof is the same as [11, Lemma 6.3] and/or [9, Theorem 3.1], so we omit it.
Proof of Theorem 6.3: By Proposition 6.2, we know that t}F pwk q}u is monotonically nonincreasing,hence convergent. Assume that t}F pwk q}u does not converge to zero. Then, there exists ϵ ą 0 suchthat twk u Ă Ωpϵq for inﬁnitely many k. We will show that the following two cases both lead tocontradictions and thus, the hypothesis }F pwk q} Û 0 is not valid.
Case 1. For inﬁnitely many k, if step (3b) in Algorithm 5 is executed with αk ” ᾱk , it follows
from Proposition 6.2 that φpwk`1 q{φpwk q ď λk :“ r1 ´ 2ᾱk β p1 ´ σk qs . Since tᾱk u is bounded awayfrom zero by Lemma 7.6 and tσk u is bounded away from one, then tλk u is bounded away from one
and hence, φpwk q Ñ 0; this is a contradiction.
Case 2. On the other hand, for inﬁnitely many k, if αk ă ᾱk , we have that αk ď θᾱk .
Then, condition (33) fails to hold for an α̃k with αk ă α̃k ď αk {θ “ θt´1 ᾱk . Notice that αk {θ isthe value corresponding
to the˘ last failure. Recall that the derivative of the real-valued function`¯ w pα∆wk q isα ÞÑ φpαq :“ φ Expk¯ w pα∆wk qqrDExp¯ w pα∆wk q r∆wk ss.
φ1 pαq “ DφpExp
k
k

18

(37)

Applying the mean value theorem to φpαq on the interval r0, α̃k s yields a number ξ P p0, 1q suchthat α̃k φ1 pξ α̃k q “ φpα̃k q ´ φp0q. For short, let u :“ ξ α̃k ∆wk . Hence,α̃k β xgrad φk , ∆wk y
ă φpα̃k q ´ φp0q (as condition (33) fails for α̃k )“ α̃k φ1 pξ α̃k q“‰¯ w puqq DExp
¯ w puq r∆wk s (by equation (37))“ α̃k DφpExp
k

(38)

k

¯ w puqq, DExp
¯ w puq r∆wk sy.
“ α̃k xgrad φpExpkk
On the other hand, note that
xgrad φk , ∆wk y“ xgrad φk , uy{ξ α̃k(39)¯ w puq rgrad φk s , DExp
¯ w puq rusy{ξ α̃k (by Lemma 4.12)“ xDExpkk¯¯ w puq r∆wk sy.
“ xDExpw puq rgrad φk s , DExpkk

Subtracting α̃k xgrad φk , ∆wk y from both sides of (38) and using equalities (39) givesα̃k pβ ´ 1q xgrad φk , ∆wk y“‰¯ w puqq, DExp
¯ w puq r∆wk sy ´ xgrad φk , ∆wk yă α̃k xgrad φpExpkk¯ w puqq ´ DExp
¯ w puq rgrad φk s , DExp
¯ w puq r∆wk sy“ α̃k xgrad φpExp
k
k
k
›
››
›
¯
¯
¯
›
›
›
ď α̃k grad φpExpwk puqq ´ DExpwk puq rgrad φk s DExpwk puq r∆wk s›››››¯ w puq rgrad φpwk qs› ›DExp¯ w puq r∆wk s›“ α̃k ›grad φpyq ´ DExpkk
¯ w puq)(by letting y :“ Expkď α̃k κ }u} }∆wk } (by (22), Deﬁnition 4.2 and Lemma 4.12)
2“ κξ α̃k2 }∆wk } .
2Finally, we obtain pβ ´ 1q xgrad φk , ∆wk y {pκξ }∆wk } q ă α̃k .
Because xgrad φk , ∆wk y ă 0 and αk satisﬁes (33), we have
φk p0q ´ φk pαk q ě ´αk β xgrad φk , ∆wk yě ´θβ α̃k xgrad φk , ∆wk y2

ě ´θβ xgrad φk , ∆wk y pβ ´ 1q xgrad φk , ∆wk y {pκξ }∆wk } q2

ě rθβp1 ´ βq{κξs pxgrad φk , ∆wk y { }∆wk }q“ ω pxgrad φk , ∆wk y { }∆wk }q ,where ωp¨q is an F -function (see [28, Deﬁnition 14.2.1 & 14.2.2 in P479]). Since tφk u is bounded below and φk ě φk`1 , it follows that limkÑ8 pφk ´ φk`1 q “ 0. By the deﬁnition of F -functions,
we obtain xgrad φk , ∆wk y { }∆wk } Ñ 0. Since t}∆wk }u is bounded ?(Proposition 7.5), we havexgrad φk , ∆wk y Ñ 0. Choosing ρk with zkT sk {m ď ρk ď }F pwk q} { m in Algorithm 5 implies
that?
xgrad φk , ∆wk y {p´2q “ φk ´ σk ρk zkT sk ě φk ´ σk }F pwk q} zkT sk { m2

ě φk ´ σk }F pwk q} ě p1 ´ σk qφk .
This shows that φpwk q Ñ 0, because tσk u is bounded away from one; this is a contradiction.

8

Numerical Experiments

The numerical experiments compared the performance of the globally convergent RIPM (Algorithm5) with those of other Riemannian methods in solving two problems. They were conducted in
Matlab R2022a on a computer with an Intel Core i7-10700 (2.90GHz) and 16GB RAM. Algorithm
5 utilized Manopt 7.0 [6], a Riemannian optimization toolbox for Matlab. The problems involve
three manifolds:19• ﬁxed-rank manifold, Mr “ tX P Rmˆn : rankpXq “ ru;␣(• Stiefel manifold, Stpn, kq “ X P Rnˆk : X T X “ Ik ;`˘• oblique manifold, Obpn, kq “ tX P Rnˆk : X T X ii “ 1, @i “ 1, . . . , ku.
We only consider their embedded geometry and we apply the default retractions in Manopt, e.g.,
the retraction based on QR decomposition for the Stiefel manifold.
Problem I. Recently, [31] proposed the Nonnegative Low-Rank Matrix (NLRM) approximation.
Formally, NLRM aims to solvemin }A ´ X}2F

XPMr

s.t. X ě 0.

(NLRM)

Input. We tested three cases of integer triples pm, n, rq: p20, 16, 2q, p30, 24, 3q, and p40, 32, 4q. For
each pm, n, rq, we generated nonnegative L P Rmˆr , R P Rrˆn whose entries follow a uniform
distribution in [0,1]. Then, we added the Gaussian noise with zero mean and diﬀerent standard
deviation pσ “ 0, 0.001, 0.01q to A “ LR. When there is no noise (i.e., σ “ 0), the input data
matrix A itself is exactly a solution.
Problem II. Given C P Rnˆk , [19] computed its projection onto the nonnegative part of the
Stiefel manifold. If the distance is measured in terms of }C ´ X}2F , we can express it equivalently as
min
XPStpn,kq´2 tracepX T Cq

s.t. X ě 0.

(Model St)

In [19], it is shown that (Model St) can be equivalently reformulated intomin
XPObpn,kq´2 tracepX T Cq

s.t. X ě 0, }XV }F “ 1,(Model Ob)

where the positive integer p and V P Rkˆp can be arbitrary as long as }V }F “ 1 and V V T is entrywise positive. We examined both models. Input. We considered four cases of integer pairs pn, kq:
p40, 8q, p50, 10q, p60, 16q, and p70, 14q. For a general C, it is always diﬃcult to seek nonnegative projections globally. Fortunately, Proposition 1 in [19] showed a way to construct C such that Model St hasa unique, known solution X ˚ . In our experiments, we generated a feasible point B of Model St; then,we set C by codes: X1=(B>0).*(1+rand(n,k));Xstar=X1./sqrt(sum(X1.*X1));L=rand(k,k);
L=L+k*eye(k);C=Xstar*L’. The initial point was computed by projecting C onto the Stiefel manifold. In addition, for Model Ob we set p=1;V=ones(k,p);V=V/norm(V,"fro").

8.1

Implementation Details

The experimental implementation of Algorithm 5 (i.e.,RIPM) initialized z0 and s0 from a uniform distribution in [0,1] and set y0 “ 0 if equality constraints exist. We used ρk “ zkT sk {m,σk “ mint0.5, }F pwk q}1{2 u and Algorithm 3 to solve the condensed form of Newton equation (19).
Algorithm 3 stopped when the relative residual went below 10´9 , or it reached 1000 iterations. We
used a backtracking line search simultaneously for the central conditions and suﬃcient decreasing
conditions. We set γ´1 “ 0.9, γk`1 “ pγk ` 0.5q{2; and β “ 10´4 , θ “ 0.5. We compared RIPM
with the following Riemannian methods:• RALM: Riemannian augmented Lagrangian method [23].
• REPMlqh : Riemannian exact penalty method with smoothing functions of linear-quadratic
and pseudo-Huber [23].
• REPMlse : Riemannian exact penalty method with smoothing functions of log-sum-exp [23].
• RSQP: Riemannian sequential quadratic programming [27].
Let rts` :“ maxp0, tq and rts´ :“ minp0, tq for any t P R. The experimental settings followed thoseof [27], where they used the KKT residual deﬁned as
g
f
m
l
f
ÿ
ÿ
2
e}grad Lpwq}2 `2 ` rg pxqs2 ` |z g pxq|2 u `trzs|hi pxq|
i ´
i
i i
`
x
i“1

i“1

20

to measure the deviation of an iterate from the KKT conditions. For the parameters of RALM,
REPMs, and RSQP, we utilized the experimental setting and Matlab codes provided by [27].
We conducted 20 random trials of each problem and model. All the algorithms ran with thesame initial point. The experiment is considered successfully terminated if it ﬁnds a solution witha KKT residual lower than ϵkkt before triggering any of the stopping conditions. For the ﬁrst-orderalgorithms (including RALM and the REPMs), the stopping conditions are: elapsed time exceedingtmax seconds, number of outer iterations exceeding 1,000, or failure of the algorithm to update any
parameters. For the second-order algorithms (including RSQP and RIPM), the stopping conditions
are elapsed time exceeding tmax seconds or a number of outer iterations exceeding 10,000. Here,considering that some problems might not have converged easily, the maximum number of iterations
was chosen to be 1,000 (10,000), which was a suﬃciently large value. The selection of tmax is related
to the actual time it took to run all the codes on the computer. Setting tmax too large resulted
in excessive time spent on poorly performing algorithms. On the other hand, ϵkkt was chosen tobetter demonstrate that second-order algorithms could achieve more accurate solutions. Therefore,
we chose the appropriate values for tmax and ϵkkt according to the problem that was to be solved.

8.2

Results and Analysis

The tables in this subsection report the success rate (Success), the average time (Time), and theaverage iteration number (Iter.) among the successful trials. In order to capture the combinationof stability and speed, the algorithms with 0.9 or higher are ﬁrst highlighted in bold in the “Success”
column, and then, of those algorithms, the fastest result is highlighted in bold in the “Time” column.
Here, “successful convergence in numerical experiments” means that the algorithm can generate a
relatively accurate solution in a relatively reasonable amount of time. It is not exactly the sameas the “theoretical convergence in any global convergence theorem”. Numerical experiments reﬂect
the actual performance of the algorithm in the application. For example, the ﬁrst-order algorithms(RALM, REPMs) also have theoretical global convergence, but under our high criterion, it is diﬃcult
for them to generate a high-precision solution within a certain period of time. The second-orderalgorithm (RIPM, RSQP), on the other hand, is excellent in terms of accuracy, although it takes
quite a bit of time.
Problem I. Here, we set tmax “ 180, ϵkkt “ 10´8 . The numerical results in Table 1 show that
RIPM had the best performance, except for cases (30, 40, 3) and (40, 32, 4) without noise. The
time spent by RALM and the REPMs grew slowly with the problem size m and n, but their successrates dropped sharply as the noise level (standard deviation σ) intensiﬁed, eventually leading to
non-convergence. In contrast, RSQP and RIPM were more stable and robust, while RIPM wasmuch faster than RSQP. The cost of RSQP increased drastically with the problem size because aquadratic programming subproblem (deﬁned over tangent space of current point) had to be solved ineach iteration. Unlike RIPM, which uses Krylov subspace methods introduced in Subsection 3.4.2 toavoid expensive computations, RSQP had to transform the subproblem into a matrix representation(similar to Step 1-6 in Subsection 3.4.2). As can be seen from Table 1, RIPM took about the sameamount of time as RALM and the REPMs did.
Problem II. Here, we set tmax “ 600, ϵkkt “ 10´6 for both (Model St) and (Model Ob). Since
the true solution X ˚ is known, we added a column showing the average error }X̃ ´ X ˚ }F , whereX̃ denotes the ﬁnal iterate. The numerical results are listed in Table 2 and 3. The Error columns
show that if the KKT residual is suﬃciently small, then X̃ does approximate the true solution. In
particular, RSQP and RIPM yield a more accurate solution (the error is less than 10 ´7 ). From
Table 2, we can see that RALM is stable and fast for (Model St). However, from Table 3, thesuccess rate of RALM for (Model Ob) decreases as the problem size becomes larger. The REPMs
do not work at all on either model. RSQP also does not perform well on either model. In contrast,
RIPM successfully solved all instances of both models and was only slightly slower than RALM insome cases. Overall, our RIPM was relatively fast and most stable.

9Conclusions

In this paper, we proposed a Riemannian version of the classical interior point method and established its local and global convergence. To our knowledge, this is the ﬁrst study to apply the
primal-dual interior point method to the nonconvex constrained optimization problem on a Riemannian manifold. Numerical experiments showed the stability and eﬃciency of our method.
21

Table 1: Performance of various Riemannian methods on problem I (NLRM).
pm, n, rq(20,16,2)(30,24,3)(40,32,4)
no noise
Success Time [s]Iter. Success Time [s]Iter. Success Time [s]
RALM
0.4
1.115
31
0.65
1.813
31
0.75
2.800
REPMlqh 1
5.165ˆ10´1
31
1
1.009
31
1
1.747
REPMlse
1
2.242
31
1
4.041
31
0.95
6.952
RSQP
0.9
6.429
7
0.9
3.944ˆ10 8
0.9
1.254ˆ102´1
1
4.920 ˆ 10
19
1
2.247
27
1
5.277
RIPM
pm, n, rq(20,16,2)(30,24,3)(40,32,4)
σ “ 0.001 Success Time [s]Iter. Success Time [s]Iter. Success Time [s]
RALM
0.2
1.001
31
0.15
2.050
31
0.05
2.758
4.983ˆ10´1
32
0.25
1.035
31
0.15
1.787
REPMlqh 0.1
REPMlse
0.15
2.444
31
0.1
4.867
31
0.05
8.371
RSQP
0.95
6.619
7
0.95
3.848ˆ10 8
0.9
1.299ˆ102´1
RIPM
1
5.376 ˆ 10
20
1
2.342
27
1
4.631
pm, n, rq(20,16,2)(30,24,3)(40,32,4)
σ “ 0.01
Success Time [s]Iter. Success Time [s]Iter. Success Time [s]
RALM
0
0
0
REPMlqh 0000
00
REPMlse
RSQP
1
7.295
8
0.95
4.114ˆ10 8
0.95
1.430ˆ102´1
RIPM
1
5.980ˆ10
21
0.95
1.883
25
0.95
4.602

Table 2: Performance of various Riemannian methods on problem II (Model St).
(40,8)(50,10)
Success Time [s]Iter. Error
Success Time [s]Iter. Error
RALM
1
2.347
45
5.41ˆ10´7 1
4.344
545.21ˆ10´7
REPMlqh 0
0
REPMlse 0
0
RSQP
0.9
1.352ˆ10 7
2.05ˆ10´9 0.7
3.097ˆ106
2.47ˆ10´9
´8
1
2.225
31
3.72ˆ101
3.785
32
3.38ˆ10´8RIPM
pn, kq(60,12)(70,14)
Success Time [s]Iter. Error
Success Time [s]Iter. Error
RALM
1
4.097
34
4.93ˆ10´7 1
6.234
37
5.34ˆ10´7
REPMlqh 0
00
REPMlse 0
RSQP
0.65
7.802ˆ10 7
6.48ˆ10´9 0.85
1.661ˆ102 7
2.64ˆ10´9
RIPM
1
5.555
32
2.81ˆ10´8 1
7.574
332.45ˆ10´8pn, kq

Table 3: Performance of various Riemannian methods on problem II (Model Ob).
(40,8)(50,10)
Success Time [s]Iter. Error
Success Time [s]Iter. Error
RALM
1
2.510
51
5.04ˆ10´7
0.95
4.727
64
4.94ˆ10´7
REPMlqh 0
00
REPMlse 0
RSQP
0.65
8.618
5
2.30ˆ10´10 0.7
2.782ˆ10 6
1.12ˆ10´10
1
3.791
22
5.62ˆ10´9
1
5.880
23
7.93ˆ10´9RIPM
pn, kq(60,12)(70,14)
Success Time [s]Iter. Error
Success Time [s]Iter. Error
´7
RALM
0.6
5.725
49
3.82ˆ10
0.6
8.223
52
3.85ˆ10´70
REPMlqh 0
REPMlse 00
0.7
4.446ˆ10 5
1.17ˆ10´90.5
9.138ˆ10 5
1.82ˆ10´9
RSQP
RIPM
1
7.134
23
9.69ˆ10´9
1
9.268
24
1.06ˆ10´8pn, kq

22

Iter.
31
31
31
8
32
Iter.
31
31
31
8
29
Iter.
9
29

Recently, Hirai et al. [15] extended the self-concordance-based interior point methods to Riemannian manifolds. They aimed to minimize a geodesically convex (i.e., convex on manifolds) objectivef : D Ñ R deﬁned on a geodesically convex subset D Ă M. In contrast, in (RCOP), we do not require any convexity. In practice, many convex functions (in the Euclidean sense) are not geodesicallyconvex on some interested manifolds. For example, for any geodesically convex function deﬁned on aconnected, compact Riemannian manifold (e.g., Stiefel manifold), it must be constant [5, Corollary
11.10], which is not of interest in the ﬁeld of optimization. Thus, (RCOP) has a wider applicability.
In closing, let us make a comparison with the Euclidean Interior Point Method (EIPM) to
illustrate the theoretical advantages of our RIPM and discuss two future directions of research on
more advanced RIPM methods.
Comparison: Riemannian IPM (RIPM) v.s. Euclidean IPM (EIPM).
1. RIPM generalizes EIPM from Euclidean space to general Riemannian manifolds. EIPM is a
special case of RIPM when M “ Rn or Rmˆn in (RCOP).
2. RIPM inherits all the advantages of Riemannian optimization. For example, we can exploit thegeometric structure of M, which is usually regarded as a set of constraints from the Euclidean
viewpoint.
3. Note that in both RIPM and EIPM, we have to solve the condensed Newton equation (19) ateach iteration. However, if the equality constraints can be considered to be a manifold, RIPM
can solve (19) with a smaller order on Tx M ˆ Rl . For example, the problem II (Model St)
can be rewritten asmin ´2 tracepX T Cq

XPRnˆks.t. X T X “ Ik , X ě 0.

Here, Stiefel manifold is replaced by the equality constraints, i.e., we deﬁne h : Rnˆk Ñ
Sympkq : X ÞÑ hpXq :“ X T X ´ Ik ; and M “ Rnˆk , l “ dim Sympkq “ kpk ` 1q{2 in(RCOP). Then, when we apply EIPM, it requires us to solve (19) of order nk ` kpk ` 1q{2.
On the other hand, if we apply RIPM to (Model St), then (19) reduces to (21) since there
are only inequality constraints on M “ Stpn, kq. In this case, we solve the equation of ordernk ´ kpk ` 1q{2, i.e., the dimension of Stpn, kq. Compared to EIPM, using RIPM reduces ourdimensionality by kpk ` 1q.
4. RIPM can solve some problems that EIPM cannot. For example, the problem I (NLRM) canbe rewritten asmin }A ´ X}2F s.t. rankpXq “ r, X ě 0.
XPRmˆnSince the rank function, X ÞÑ rankpXq, is not even continuous, we cannot apply EIPM.
Future Work I: Preconditioner for linear operator equation. With regard to the complementary condition, Sk´1 Zk values display a huge diﬀerence in magnitude as k Ñ 8. The operator
Θ :“ Gx S ´1 ZG˚x causes the system (19) to be ill-conditioned, risking failure of the iterative method
without preconditioning. Matrix-decomposition-based preconditioner methods cannot be applied toa problem that does not have a matrix form. A possible alternative is to ﬁnd a nonsingular P suchthat the condition number of P ´1 T is smaller.
Future Work II: Treatment of more state-of-the-art interior point methods. While
we have considered interior point methods on a manifold for the ﬁrst time, our Euclidean theoretic
counterpart is an early nonlinear interior point method algorithm [11]; however, the counterpart now
appears to be obsolete compared with more recent interior point methods. For example, our method
does not drive the iteration towards minimizers but only towards stationary points; globalization isdone by monitoring only the KKT residuals; moreover, the boundedness assumption (C2) of tzk u
is too strong to hold in some simple cases (see Wächter-Biegler eﬀect [32]). It remains an important
issue to adapt more modern interior point methods to manifolds, although we may encounter various
diﬃculties in Riemannian geometry.
Acknowledgements We are deeply grateful to Hiroyuki Sato of Kyoto University for his valuable comments and suggestions. This work was supported by JST SPRING, Grant Number JPMJSP2124, and JSPS KAKENHI, Grant Numbers 19H02373 and 22K18866.
23

Data availability The datasets generated during and/or analyzed during the current study areavailable from the corresponding author upon reasonable request.

References[1] P.-A. Absil, R. Mahony, and R. Sepulchre, Optimization Algorithms on Matrix Manifolds,
Princeton University Press, Princeton, NJ, 2009.
[2] R. Bergmann and R. Herzog, Intrinsic formulation of KKT conditions and constraint qualiﬁcations on smooth manifolds, SIAM J. Optim., 29 (2019), pp. 2423–2444.
[3] S. Bonettini, E. Galligani, and V. Ruggiero, An inexact Newton method combined with
Hestenes multipliers scheme for the solution of Karush–Kuhn–Tucker systems, Appl. Math.
Comput., 168 (2005), pp. 651–676.
[4] M.A.d.A. Bortoloti, T.A. Fernandes, O.P. Ferreira, and J. Yuan, Damped Newton s method
on Riemannian manifolds, J. Glob. Optim., 77 (2020), pp. 643–660.
[5] N. Boumal, An Introduction to Optimization on Smooth Manifolds, Cambridge University
Press, (2023).
[6] N. Boumal, B. Mishra, P.-A. Absil, and R. Sepulchre, Manopt, a Matlab toolbox for optimization on manifolds, J. Mach. Learn. Res., 15 (2014), pp. 1455–1459.
[7] M.P. Do Carmo, Riemannian Geometry, Birkhäuser, Boston, MA, 1992.
[8] J.E. Dennis, R.B. Schnabel, Numerical Methods for Unconstrained Optimization and Nonlinear Equations, SIAM, 1996.
[9] C. Durazzi, On the Newton interior-point method for nonlinear programming problems, J.
Optim. Theory Appl., 104 (2000), pp. 73–90.
[10] C. Durazzi and V. Ruggiero, Global convergence of the Newton interior-point method for
nonlinear programming, J. Optim. Theory Appl., 120 (2004), pp. 199–208.
[11] A. El-Bakry, R.A. Tapia, T. Tsuchiya, and Y. Zhang, On the formulation and theory of
the Newton interior-point method for nonlinear programming, J. Optim. Theory Appl., 89(1996), pp. 507–541.
[12] T.A. Fernandes, O.P. Ferreira, and J. Yuan, On the superlinear convergence of Newton s
method on Riemannian manifolds, J. Optim. Theory Appl., 173 (2017), pp. 828–843.
[13] O.P. Ferreira and B.F. Svaiter, Kantorovich’s theorem on Newton’s method in Riemannian
manifolds, J. Complex., 18 (2002), pp. 304–329.
[14] K.A. Gallivan, C. Qi, and P.-A. Absil, A Riemannian Dennis-Moré condition, in HighPerformance Scientiﬁc Computing: Algorithms and Applications, M.W. Berry, K.A. Gallivan, E. Gallopoulos, A. Grama, B. Philippe, Y. Saad and F. Saied, eds., Springer, London,
2012, pp. 281–293.
[15] H. Hirai, H. Nieuwboer, and M. Walter, Interior-point methods on manifolds: theory andapplications, arXiv:2303.04771, (2023) to appear in FOCS 2023.
[16] J. Hu, X. Liu, Z.W. Wen, and Y.X. Yuan, A brief introduction to manifold optimization,
J. Oper. Res. Soc. China, 8 (2020), pp. 199–248.
[17] W. Huang, Optimization Algorithms on Riemannian Manifolds with Applications, Ph.D.
thesis, Florida State University, Tallahassee, FL, 2013.
[18] W. Huang, P.-A. Absil, and K.A. Gallivan, A Riemannian symmetric rank-one trust-region
method, Math. Program., 150 (2015), pp. 179–216.
[19] B. Jiang, X. Meng, Z.W. Wen, and X. Chen, An exact penalty approach for optimization
with nonnegative orthogonality constraints, Math. Program., (2022), pp. 1–43.
24[20] M. Kojima, S. Mizuno, and A. Yoshise, A primal-dual interior point algorithm for linear
programming, in Progress in Mathematical Programming, N. Megiddo ed., Springer, 1989,
pp. 29–47.
[21] J.M. Lee, Introduction to Smooth Manifolds, 2nd ed., Grad. Texts in Math. 218, Springer,
2013.
[22] J.M. Lee, Introduction to Riemannian Manifolds, 2nd ed., Grad. Texts in Math. 176,
Springer, 2018.
[23] C. Liu and N. Boumal, Simple algorithms for optimization on Riemannian manifolds withconstraints, Appl. Math. Optim., 82 (2020), pp. 949–981.
[24] I. Lustig, R.E. Marsten, and D.F. Shanno, Computational experience with a primal-dual
interior point method for linear programming, Linear Algebra Appl., 152 (1991), pp. 191–222.
[25] N. Megiddo, Pathways to the optimal set in linear programming, in Progress in Mathematical
Programming, N. Megiddo ed., Springer, 1989, pp. 131–158.
[26] J. Nocedal and S.J. Wright, 2nd ed., Numerical Optimization, Springer, New York, NY,
2006.
[27] M. Obara, T. Okuno, and A. Takeda, Sequential quadratic optimization for nonlinear optimization problems on Riemannian manifolds, SIAM J. Optim., 32 (2022), pp. 822–853.
[28] J.M. Ortega and W.C. Rheinboldt, Iterative Solution of Nonlinear Equations in Several
Variables, SIAM, 2000.
[29] Y. Saad, 2nd ed., Iterative Methods for Sparse Linear Systems, Other Titles in Applied
Mathematics SIAM, 2003.
[30] A. Schiela and J. Ortiz, An SQP method for equality constrained optimization on Hilbert
manifolds, SIAM J. Optim., 31 (2021), pp. 2255–2284.
[31] G.J. Song and M.K. Ng, Nonnegative low rank matrix approximation for nonnegative matrices, Appl. Math. Lett., 105 (2020), p. 106300.
[32] A. Wächter and L.T. Biegler, Failure of global convergence for a class of interior point
methods for nonlinear programming, Math. Program., 88 (2000), pp. 565–574.
[33] S.J. Wright, Primal-dual Interior-point Methods, SIAM, 1997.
[34] Y. Yamakawa and H. Sato, Sequential optimality conditions for nonlinear optimization on
Riemannian manifolds and a globally convergent augmented Lagrangian method, Comput.
Optim. Appl., 81 (2022), pp. 397–421.
[35] H. Yamashita and H. Yabe, Superlinear and quadratic convergence of some primal-dual
interior point methods for constrained optimization, Math. Program., 75 (1996), pp. 377–
397.
[36] W.H. Yang, L.H. Zhang, and R. Song, Optimality conditions for the nonlinear programming
problems on Riemannian manifolds, Paciﬁc J. Optim., 10 (2014), pp. 415–434.
[37] Y. Ye, Interior Point Algorithms: Theory and Analysis, John Wiley & Sons, 1997.
[38] X. Zhu and H. Sato, Riemannian conjugate gradient methods with inverse retraction, Comput. Optim. Appl., 77 (2020), pp. 779–810.

25