# Riemannian Interior Point Methods
for Constrained Optimization on Manifolds111An earlier version of this article
has been circulated under the title “Superlinear and Quadratic Convergence of
Riemannian Interior Point Methods.”

Zhijian Lai  Akiko Yoshise  Zhijian Lai  University of Tsukuba, Tsukuba,
Ibaraki 305-8573, Japan<EMAIL_ADDRESS>Akiko Yoshise  Corresponding
author. University of Tsukuba, Tsukuba, Ibaraki 305-8573, Japan.
<EMAIL_ADDRESS>
(Revised )

###### 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. Numerical experiments show the stability and efficiency of our
method.

Communicated by Sándor Zoltán Németh

Keywords: Riemannian manifolds, Riemannian optimization, Nonlinear
optimization and Interior point method.

##  1 Introduction

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

$\begin{array}[]{ll}\underset{x\in\mathcal{M}}{\min}&f(x)\\\ \text{ s.t.
}&h(x)=0,\text{ and }g(x)\leq 0,\\\ \end{array}$ (RCOP)

where $\mathcal{M}$ is a connected, complete $d$-dimensional Riemannian
manifold and
$f\colon\mathcal{M}\to\mathbb{R},h\colon\mathcal{M}\to\mathbb{R}^{l}$, and
$g\colon\mathcal{M}\to\mathbb{R}^{m}$ are smooth functions. This problem
appears in many applications, for instance, matrix approximation with
nonnegative constraints on a fixed-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 back to
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 Qualifications (CQs) on manifolds. Yamakawa and Sato [34] proposed
sequential optimality conditions in the Riemannian case. Liu and Boumal [23]
were the first to develop practical algorithms. They extended the augmented
Lagrangian method and exact penalty method to (RCOP). Schiela and Ortiz [30]
and Obara 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 field
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 first decade of the 21st
century, a large number of interior point methods for nonlinear programming
emerged. They proved to be as successful as the linear ones. A subclass known
as _primal-dual_ interior point methods is the most efficient and practical.
As described in [24], the primal-dual approach to linear programming was
introduced in [25]: it was first 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., $\mathcal{M}=\mathbb{R}^{d}$ in
(RCOP), to the Riemannian setting. We call this extension the _Riemannian
Interior Point Method_ (RIPM). Our contributions are summarized as follows:

  1. 1.

To our knowledge, this is the first study to apply the primal-dual interior
point method to the nonconvex constrained optimization problem on Riemannian
manifolds. One significant contribution is that we establish many essential
foundational concepts for the general interior point method in the Riemannian
context, such as the KKT vector field and its covariant derivative. In
addition, we build the first framework for the Riemannian version of the
interior point method. These contributions will have uses in the future,
especially in developing more advanced interior point methods.

  2. 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 for nonlinear nonconvex
programming first proposed by El-Bakry et al. [11].

  3. 3.

Our numerical experiments222The code is freely available at
https://doi.org/10.5281/zenodo.10612799 demonstrate the great potential of
RIPM. The method meets the challenges presented in these experiments with
better stability and higher accuracy compared with 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 the notation of Riemannian geometry and explain the Riemannian Newton
method. In Section 3, we give a full interpretation of our RIPM and describe a
prototype algorithm of RIPM. We also investigate the use of Krylov subspace
methods to efficiently 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 a
classical line search; then, Section 7 proves its global convergence. Section
8 is a collection of numerical experiments. Section 9 summarizes our research
and presents future work.

##  2 Notation

###  2.1 Riemannian Geometry

Let us briefly review some concepts from Riemannian geometry, following the
notation of [5]. $\mathcal{M}$ denotes a finite-dimensional smooth manifold.
Let $p\in\mathcal{M}$ and $T_{p}\mathcal{M}$ be the tangent space at $p$ with
$0_{p}$ being its zero element. We use a canonical identification
$T_{p}\mathcal{E}\cong\mathcal{E}$ for a vector space $\mathcal{E}$ and
$p\in\mathcal{E}$. A vector field is a map $V\colon\mathcal{M}\to
T\mathcal{M}$ with $V(p)\in T_{p}\mathcal{M}$, where
$T\mathcal{M}:=\bigcup_{p\in\mathcal{M}}T_{p}\mathcal{M}$ is the tangent
bundle. $\mathfrak{X}(\mathcal{M})$ denotes the set of all smooth vector
fields defined on $\mathcal{M}$. Furthermore, $\mathcal{M}$ is a Riemannian
manifold if it is equipped with a Riemannian metric, that is, a choice of
inner product $\langle\cdot,\cdot\rangle_{p}\colon T_{p}\mathcal{M}\times
T_{p}\mathcal{M}\to\mathbb{R}$ for each tangent space at $p$ on $\mathcal{M}$
such that for all $V,W\in\mathfrak{X}(\mathcal{M})$, the map

$p\mapsto\langle V(p),W(p)\rangle_{p}$ (1)

is a smooth function from $\mathcal{M}$ to $\mathbb{R}$. Riemannian metric
induces the norm $\|\xi\|_{p}:=\sqrt{\langle\xi,\xi\rangle_{p}}$ for $\xi\in
T_{p}\mathcal{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 $\mathcal{M}$,
$c:[a,b]\rightarrow\mathcal{M}$, the length of $c$ is defined as
$L(c):=\int_{a}^{b}\|\dot{c}(t)\|_{c(t)}\mathrm{d}t,$ where $\dot{c}(t)\in
T_{c(t)}\mathcal{M}$ is the velocity vector of $c$ at $t$. Since $\mathcal{M}$
is connected, there exists a curve segment connecting any pair of points
$p,q\in\mathcal{M}$. Indeed, $\mathcal{M}$ is a metric space under the
Riemannian distance $d(p,q):=\inf_{c}L(c)$ where the infimum is taken over all
curve segments joining $p$ to $q$ [22, Thm. 2.55]. For two manifolds
$\mathcal{M}_{1},\mathcal{M}_{2}$ and a smooth map
$F\colon\mathcal{M}_{1}\to\mathcal{M}_{2}$, the differential of $F$ at
$p\in\mathcal{M}_{1}$ is a linear operator denoted as $\mathcal{D}F(p)\colon
T_{p}\mathcal{M}_{1}\to T_{F(p)}\mathcal{M}_{2}$. Let
$\mathfrak{F}(\mathcal{M})$ be the set of all smooth _scalar fields_ (or say
real-valued functions) $f\colon\mathcal{M}\to\mathbb{R}$. The Riemannian
gradient of $f$ at $p$, $\operatorname{grad}f(p)$, is defined as the unique
element of $T_{p}\mathcal{M}$ that satisfies
$\langle\xi,\operatorname{grad}f(p)\rangle_{p}=\mathcal{D}f(p)[\xi]$ for all
$\xi\in T_{p}\mathcal{M}$, where $\mathcal{D}f(p)\colon T_{p}\mathcal{M}\to
T_{f(p)}\mathbb{R}\cong\mathbb{R}$. Note that for any
$f\in\mathfrak{F}(\mathcal{M})$, the gradient vector field
$x\mapsto\operatorname{grad}f(x)$ is a smooth vector field on $\mathcal{M}$,
i.e., $\operatorname{grad}f\in\mathfrak{X}(\mathcal{M})$. A retraction
$\mathrm{R}\colon T\mathcal{M}\to\mathcal{M}$ is a smooth map such that
$\mathrm{R}_{p}\left(0_{p}\right)=p$ and
$\mathcal{D}\mathrm{R}_{p}\left(0_{p}\right)=\operatorname{id}_{T_{p}\mathcal{M}}$,
i.e., the identity map on $T_{p}\mathcal{M}$, where $\mathrm{R}_{p}$ is the
restriction of $\mathrm{R}$ to $T_{p}\mathcal{M}$ and
$\mathcal{D}\mathrm{R}_{p}\left(0_{p}\right)\colon
T_{0_{p}}(T_{p}\mathcal{M})\cong T_{p}\mathcal{M}\to T_{p}\mathcal{M}$ is the
differential of $\mathrm{R}_{p}$ at $0_{p}$. One theoretically perfect type of
retraction is the exponential map, denoted as $\operatorname{Exp}$. Since
$\mathcal{M}$ is complete, the exponential map is well-defined on the whole
tangent bundle. Let $\operatorname{Exp}_{p}\colon
T_{p}\mathcal{M}\rightarrow\mathcal{M}$ be the exponential map at $p$; then
$t\mapsto\operatorname{Exp}_{p}(t\xi)$ is the unique geodesic that passes
through $p$ with velocity $\xi\in T_{p}\mathcal{M}$ when $t=0$.

###  2.2 Riemannian Newton Method

The Newton method is a powerful tool for finding the zeros of nonlinear
functions in the Euclidean setting. The generalized Newton method has been
studied in the Riemannian setting; it aims to find a _singularity_ for the
vector field $F\in\mathfrak{X}(\mathcal{M})$, i.e., a point $p\in\mathcal{M}$
such that,

$F(p)=0_{p}\in T_{p}\mathcal{M}.$ (2)

Let $\nabla$ be the Levi-Civita connection on $\mathcal{M}$, i.e., the unique
symmetric connection compatible with the Riemannian metric. The covariant
derivative $\nabla F$ assigns each point $p\in\mathcal{M}$ a linear operator
$\nabla F(p)\colon T_{p}\mathcal{M}\to T_{p}\mathcal{M}$. In particular, the
Riemannian Hessian of $f\in\mathfrak{F}(\mathcal{M})$ at $p$ is a self-adjoint
operator on $T_{p}\mathcal{M}$, defined as
$\operatorname{Hess}f(p):=\nabla(\operatorname{grad}f)(p).$ The Riemannian
Newton method for solving (2) is performed as Algorithm 1.

Input: A vector field $F\in\mathfrak{X}(\mathcal{M})$, an initial point
$p_{0}\in\mathcal{M}$ and a retraction $\mathrm{R}$ on $\mathcal{M}$.

Output: Sequence $\\{p_{k}\\}\subset\mathcal{M}$ such that
$\left\\{p_{k}\right\\}\rightarrow p^{*}$ and $F(p^{*})=0_{p^{*}}$.

Set $k\to 0$;

while  _Stopping criterion not satisfied_ do

1\. Obtain $\xi_{k}\in T_{p_{k}}\mathcal{M}$ by solving the Newton equation (a
linear operator equation on tangent space $T_{p_{k}}\mathcal{M}$):  $\nabla
F(p_{k})[\xi_{k}]=-F(p_{k});$ (3)

2\. Compute the next point as $p_{k+1}:=\mathrm{R}_{p_{k}}(\xi_{k})$;

3\. $k\to k+1$;

end while

Algorithm 1 Riemannian Newton Method for (2)

##  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=\operatorname{diag}(z_{1},\ldots,z_{m})$ with $z\in\mathbb{R}^{m}$. $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

$\mathcal{L}(x,y,z):=f(x)+\sum_{i=1}^{l}y_{i}h_{i}(x)+\sum_{i=1}^{m}z_{i}g_{i}(x),$

where $y\in\mathbb{R}^{l}$ and $z\in\mathbb{R}^{m}$ are Lagrange multipliers
corresponding to the equality and inequality constraints, respectively. With
respect to the variable $x$, $\mathcal{L}(\cdot,y,z)$ is a real-valued
function defined on $\mathcal{M}$, and its Riemannian gradient is

$\operatorname{grad}_{x}\mathcal{L}(x,y,z)=\operatorname{grad}f(x)+\sum_{i=1}^{l}y_{i}\operatorname{grad}h_{i}(x)+\sum_{i=1}^{m}z_{i}\operatorname{grad}g_{i}(x),$

where $\\{\operatorname{grad}h_{i}(x)\\}_{i=1}^{l}$ and
$\\{\operatorname{grad}g_{i}(x)\\}_{i=1}^{m}$ are the gradients of the
component functions of $h$ and $g$. The Riemannian KKT conditions (e.g., see
[23, Definition 2.3] or [36]) for (RCOP) are given by

$\operatorname{grad}_{x}\mathcal{L}(x,y,z)=0_{x};\;h(x)=0,g(x)\leq 0,z\geq
0;\;Zg(x)=0.$ (4)

The above conditions can be written in terms of slack variables
$\mathbb{R}^{m}\ni s:=-g(x)$, as

$F(w):=\left(\begin{array}[]{l}\operatorname{grad}_{x}\mathcal{L}(x,y,z)\\\
h(x)\\\ g(x)+s\\\ ZSe\end{array}\right)=0_{w}=\left(\begin{array}[]{l}0_{x}\\\
0\\\ 0\\\ 0\end{array}\right)\in T_{w}\mathcal{N},$ (5)

and $(z,s)\geq 0,$ where
$w:=(x,y,z,s)\in\mathcal{N}:=\mathcal{M}\times\mathbb{R}^{l}\times\mathbb{R}^{m}\times\mathbb{R}^{m}$.
Here, we have generated a vector field $F$ in (5) on the product manifold
$\mathcal{N}$, i.e., for any $w\in\mathcal{N}$, $F(w)\in T_{w}\mathcal{N}\cong
T_{x}\mathcal{M}\times\mathbb{R}^{l}\times\mathbb{R}^{m}\times\mathbb{R}^{m}.$
We will call $F\colon\mathcal{N}\to T\mathcal{N}$ in (5) above the _KKT vector
field_ for (RCOP). Note that for
$\xi=\left(\xi_{x},\xi_{y},\xi_{z},\xi_{s}\right)\in T_{w}\mathcal{N}$,

$\bar{\mathrm{R}}_{w}(\xi):=(\mathrm{R}_{x}(\xi_{x}),y+\xi_{y},z+\xi_{z},s+\xi_{s})$
(6)

is a well-defined retraction on the product manifold $\mathcal{N}$ as long as
$\mathrm{R}$ is a retraction on $\mathcal{M}$.

Now, we aim to find a singularity of the KKT vector field $F$ on $\mathcal{N}$
under some nonnegative constraints $(z,s)\geq 0$. If the Riemannian Newton
method is to be applied to (5), we must formulate the covariant derivative of
$F$ at an arbitrary $w\in\mathcal{N}$. Before that, we need some new symbols.
Fixing a point $x\in\mathcal{M}$, we can define two linear operators
$\mathcal{H}_{x}\colon\mathbb{R}^{l}\rightarrow T_{x}\mathcal{M}$ and
$\mathcal{G}_{x}\colon\mathbb{R}^{m}\rightarrow T_{x}\mathcal{M}$ by

$\mathcal{H}_{x}[v]:=\sum_{i=1}^{l}v_{i}\operatorname{grad}h_{i}(x),\quad\mathcal{G}_{x}[u]:=\sum_{i=1}^{m}u_{i}\operatorname{grad}g_{i}(x),$
(7)

respectively. Then, for $\mathcal{H}_{x}$, its adjoint operator
$\mathcal{H}_{x}^{*}\colon T_{x}\mathcal{M}\rightarrow\mathbb{R}^{l}$ is given
by

$\mathcal{H}_{x}^{*}[\xi]=\left(\left\langle\operatorname{grad}h_{1}(x),\xi\right\rangle_{x},\cdots,\left\langle\operatorname{grad}h_{l}(x),\xi\right\rangle_{x}\right)^{T}.$

Also, $\mathcal{G}_{x}^{*}\colon T_{x}\mathcal{M}\rightarrow\mathbb{R}^{m}$
with
$\mathcal{G}_{x}^{*}[\xi]=\left(\left\langle\operatorname{grad}g_{1}(x),\xi\right\rangle_{x},\cdots,\left\langle\operatorname{grad}g_{m}(x),\xi\right\rangle_{x}\right)^{T}$.

Now, using the solutions of Exercises 5.4 and 5.13 on covariant derivatives of
vector fields on product manifolds in monograph [5], the following results can
be easily derived. Given KKT vector field $F$ in (5) and $w\in\mathcal{N}$,
the covariant derivative of $F$ at $w$ is the operator $\nabla F(w)\colon
T_{w}\mathcal{N}\to T_{w}\mathcal{N}$ given by

$\nabla F(w)[\Delta
w]=\left(\begin{array}[]{l}\operatorname{Hess}_{x}\mathcal{L}(w)[\Delta
x]+\mathcal{H}_{x}[\Delta y]+\mathcal{G}_{x}[\Delta z]\\\
\mathcal{H}_{x}^{*}[\Delta x]\\\ \mathcal{G}_{x}^{*}[\Delta x]+\Delta s\\\
Z\Delta s+S\Delta z\end{array}\right),$ (8)

where $\Delta w=(\Delta x,\Delta y,\Delta z,\Delta s)\in T_{w}\mathcal{N}$ and
$\operatorname{Hess}_{x}\mathcal{L}(w)$ is the Riemannian Hessian of real-
valued function $\mathcal{L}(\cdot,y,z)$. Indeed,
$\operatorname{Hess}_{x}\mathcal{L}(w)\colon T_{x}\mathcal{M}\to
T_{x}\mathcal{M}$ satisfies

$\operatorname{Hess}_{x}\mathcal{L}(w)=\operatorname{Hess}f(x)+\sum_{i=1}^{l}y_{i}\operatorname{Hess}h_{i}(x)+\sum_{i=1}^{m}z_{i}\operatorname{Hess}g_{i}(x),$
(9)

where $\\{\operatorname{Hess}h_{i}(x)\\}_{i=1}^{l}$ and
$\\{\operatorname{Hess}g_{i}(x)\\}_{i=1}^{m}$ 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 field given in (8). Let $\mathbb{A}(x):=\left\\{i\mid
g_{i}(x)=0\right\\}$ denote the active set at $x\in\mathcal{M}$. A prior study
on Riemannian optimality conditions [36, 2] showed that the following
assumptions are meaningful. We call them the _standard Riemannian assumptions_
for (RCOP). Note that the $x^{*}$ and $w^{*}$ in (A2)-(A4) all refer to those
in (A1).

  1. (A1)

Existence. There exists $w^{*}=(x^{*},y^{*},z^{*},s^{*})$ satisfying the
Riemannian KKT conditions (4). Here, we introduce the slack variables
$s^{*}:=-g(x^{*})$.

  2. (A2)

Linear Independence Constraint Qualification (LICQ).
The set
$\left\\{\operatorname{grad}h_{i}(x^{*})\right\\}_{i=1}^{l}\bigcup\left\\{\operatorname{grad}g_{i}(x^{*})\right\\}_{i\in\mathbb{A}(x^{*})}$
is linearly independent in $T_{x^{*}}\mathcal{M}$.

  3. (A3)

Strict complementarity. $(z^{*})_{i}>0$ if $g_{i}(x^{*})=0$ for all
$i=1,\cdots,m$.

  4. (A4)

Second-order sufficiency.
$\left\langle\operatorname{Hess}_{x}\mathcal{L}(w^{*})[\xi],\xi\right\rangle>0$
for all nonzero $\xi\in T_{x^{*}}\mathcal{M}$
satisfying $\left\langle\xi,\operatorname{grad}h_{i}(x^{*})\right\rangle=0$
for $i=1,\dots,l,$ and
$\left\langle\xi,\operatorname{grad}g_{i}(x^{*})\right\rangle=0$ for
$i\in\mathbb{A}(x^{*})$.

Recall the Riemannian Newton method discussed in Section 2.2. It can be shown
that if $p^{*}$ is a solution of equation (2) and the covariant derivative
$\nabla F\left(p^{*}\right)$ is _nonsingular_ , then Algorithm 1 has the local
superlinear convergence [12] and local quadratic convergence [13] under
certain mild conditions on the map $p\mapsto\nabla F(p)$. 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^{*}$,
then the operator $\nabla F(w^{*})$ in (8) is nonsingular.

###### Proof.

This proof omits all the asterisks of the variables. Define
$\mathbb{E}:=\\{1,\dots,l\\}$ and $\mathbb{I}:=\\{1,\dots,m\\}$. Take some
$w=(x,y,z,s)\in\mathcal{N}$ satisfying (A1)-(A4), then we have
$s_{i}=-g_{i}(x)$ and $z_{i}s_{i}=0$ for all $i\in\mathbb{I}$. For short, let
$\mathbb{A}:=\mathbb{A}(x)\subset\mathbb{I}$. Suppose that
$\nabla{F}(w)[\Delta w]=0$ for some $\Delta w=(\Delta x,\Delta y,\Delta
z,\Delta s)\in T_{w}\mathcal{N}\cong
T_{x}\mathcal{M}\times\mathbb{R}^{l}\times\mathbb{R}^{m}\times\mathbb{R}^{m}$.
$\Delta y_{i}$ denotes the components of the vector $\Delta y$, as do $\Delta
z_{i}$, $\Delta s_{i}$. To prove its nonsingularity, we will show that $\Delta
w=0$. Expanding the equation $\nabla{F}(w)[\Delta w]=0$ gives

$\left\\{\begin{aligned} 0&=\operatorname{Hess}_{x}\mathcal{L}(w)[\Delta
x]+\sum_{i\in\mathbb{E}}\Delta
y_{i}\operatorname{grad}h_{i}(x)+\sum_{i\in\mathbb{I}}\Delta
z_{i}\operatorname{grad}g_{i}(x),\\\
0&=\left\langle\operatorname{grad}h_{i}(x),\Delta x\right\rangle,\text{ for
all }i\in\mathbb{E},\\\ 0&=\left\langle\operatorname{grad}g_{i}(x),\Delta
x\right\rangle+\Delta s_{i},\text{ for all }i\in\mathbb{I},\\\ 0&=z_{i}\Delta
s_{i}+s_{i}\Delta z_{i},\text{ for all }i\in\mathbb{I}.\end{aligned}\right.$
(10)

Strict complementarity (A3) and the last equalities above imply that $\Delta
s_{i}=0$ for all $i\in\mathbb{A}$ and $\Delta z_{i}=0$ for
$i\in\mathbb{I}\setminus\mathbb{A}$. Substituting those values into the system
(10) reduces it to

$\left\\{\begin{aligned} 0&=\operatorname{Hess}_{x}\mathcal{L}(w)[\Delta
x]+\sum_{i\in\mathbb{E}}\Delta
y_{i}\operatorname{grad}h_{i}(x)+\sum_{i\in\mathbb{A}}\Delta
z_{i}\operatorname{grad}g_{i}(x),\\\
0&=\left\langle\operatorname{grad}h_{i}(x),\Delta x\right\rangle,\text{ for
all }i\in\mathbb{E},\\\ 0&=\left\langle\operatorname{grad}g_{i}(x),\Delta
x\right\rangle,\text{ for all }i\in\mathbb{A},\end{aligned}\right.$ (11)

and $\Delta s_{i}=-\left\langle\operatorname{grad}g_{i}(x),\Delta
x\right\rangle\text{ for all }i\in\mathbb{I}\setminus\mathbb{A}.$ It follows
from system (11) that

$\displaystyle 0$
$\displaystyle=\langle\operatorname{Hess}_{x}\mathcal{L}(w)[\Delta
x]+\sum_{i\in\mathbb{E}}\Delta
y_{i}\operatorname{grad}h_{i}(x)+\sum_{i\in\mathbb{A}}\Delta
z_{i}\operatorname{grad}g_{i}(x),\Delta x\rangle$
$\displaystyle=\left\langle\operatorname{Hess}_{x}\mathcal{L}(w)[\Delta
x],\Delta x\right\rangle+\sum_{i\in\mathbb{E}}\Delta
y_{i}\left\langle\operatorname{grad}h_{i}(x),\Delta
x\right\rangle+\sum_{i\in\mathbb{A}}\Delta
z_{i}\left\langle\operatorname{grad}g_{i}(x),\Delta x\right\rangle$
$\displaystyle=\left\langle\operatorname{Hess}_{x}\mathcal{L}(w)[\Delta
x],\Delta x\right\rangle.$

Thus, from second-order sufficiency (A4), $\Delta x$ must be zero elements.
And then $\Delta s_{i}=0\text{ for all }i\in\mathbb{I}\setminus\mathbb{A}.$
Next, substituting $\Delta x=0$ into the first equation in (11) yields
$0=\sum_{i\in\mathbb{E}}\Delta
y_{i}\operatorname{grad}h_{i}(x)+\sum_{i\in\mathbb{A}}\Delta
z_{i}\operatorname{grad}g_{i}(x).$ The LICQ (A2) implies that the coefficients
$\Delta y_{i}$ for $i\in\mathbb{E}$ and $\Delta z_{i}$ for $i\in\mathbb{A}$
must be zero. This completes the proof. ∎

###  3.3 Prototype Algorithm of RIPM

Applying the Riemannian Newton method directly to the KKT vector field
$F:\mathcal{N}\rightarrow T\mathcal{N}$ results in the following Newton
equation (see (3) without iteration count $k$) at each iteration:

$\nabla F(w)[\Delta w]+F(w)=0.$ (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 sufficient distance from the
boundary, we introduce a perturbed complementary equation with some barrier
parameter $\mu>0$ and define the perturbed KKT vector field:

$F_{\mu}(w):=F(w)-\mu\hat{e},\text{ and
}\hat{e}:=\hat{e}(w):=\left(0_{x},0,0,e\right).$ (13)

Notice that the perturbation term $\hat{e}$, indeed, is a special vector field
on $\mathcal{N}$, not a constant, because $0_{x}$ is essentially dependent on
$w$ and/or $x$. In fact, the covariant derivative of the perturbed KKT vector
field is the same as that of the original. From the linearity of the
connection $\nabla$, for any $w\in\mathcal{N}$ and any $\mu>0$, we have

$\nabla F_{\mu}(w)=\nabla F(w)-\mu\nabla\hat{e}(w)=\nabla F(w),$ (14)

where the last equity comes from $\nabla\hat{e}(w)[\Delta
w]=\left(0_{x},0,0,0\right)$ for all $\Delta w\in T_{w}\mathcal{N}$. Applying
the Riemannian Newton method to perturbed KKT vector field $F_{\mu}(w)$ yields
the _perturbed Newton equation_ : $\nabla F_{\mu}(w)[\Delta w]+F_{\mu}(w)=0$.
From (13) and (14), this equation is equivalent to

$\nabla F(w)[\Delta w]+F(w)=\mu\hat{e},$ (15)

which reduces to the ordinary Newton equation (12) as $\mu\to 0$. At this
point, we can describe a prototype of the Riemannian Interior Point Method
(RIPM) in Algorithm 2.

Input: An initial point
$w_{0}=\left(x_{0},y_{0},z_{0},s_{0}\right)\in\mathcal{N}$ with
$\left(z_{0},s_{0}\right)>0$ and a retraction $\mathrm{R}$ on $\mathcal{M}$.
$\hat{\gamma}\in(0,1)$, $\mu_{0}>0$.

Output: Sequence $\left\\{w_{k}\right\\}\subset\mathcal{N}$ such that
$\left\\{w_{k}\right\\}\rightarrow w^{*}$ and $w^{*}$ satisfies the KKT
conditions (4).

Set $k\to 0$;

while  _Stopping criterion not satisfied_ do

1\. Obtain $\Delta w_{k}=(\Delta x_{k},\Delta y_{k},\Delta z_{k},\Delta
s_{k})\in T_{w_{k}}\mathcal{N}$ by solving the perturbed Newton equation on
$T_{w_{k}}\mathcal{N}$:  $\nabla F(w_{k})[\Delta
w_{k}]=-F(w_{k})+\mu_{k}\hat{e};$ (16)

2\. Choose $\gamma_{k}\in[\hat{\gamma},1]$ and compute the step size
$\alpha_{k}$ defined by
$\min\left\\{1,\gamma_{k}\min_{i}\left\\{-\frac{(s_{k})_{i}}{(\Delta
s_{k})_{i}}\mid(\Delta
s_{k})_{i}<0\right\\},\gamma_{k}\min_{i}\left\\{-\frac{(z_{k})_{i}}{(\Delta
z_{k})_{i}}\mid(\Delta z_{k})_{i}<0\right\\}\right\\};$ (17)

3\. Compute the next point as
$w_{k+1}:=\bar{\mathrm{R}}_{w_{k}}\left(\alpha_{k}\Delta w_{k}\right)$, see
(6);

4\. Choose $0<\mu_{k+1}<\mu_{k}$;

5\. $k\to k+1$;

end while

Algorithm 2 Prototype Algorithm of RIPM for (RCOP)

In the step 2 of Algorithm 2, we just wish to compute a step size
$0<\alpha_{k}\leq 1$ to ensure new point $w_{k+1}$ with
$\left(z_{k+1},s_{k+1}\right)>0$. There are many schemes to achieve this
purpose. The scheme in (17) is simple but sufficient to guarantee the local
convergence, as will be proved in Section 5.

###  3.4 Solving Perturbed Newton Equation Efficiently

The challenge of Algorithm 2 is how to solve the Newton equation (16) in an
efficient manner. In this subsection, we will do this in two steps: the first
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 effort 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\in\mathcal{N}$ with $(z,s)>0$, for the KKT vector field
$F(w)$ in (5), we denote its components by $F_{x},F_{y},F_{z},F_{s}$ in top-
to-bottom order, namely,
$F_{x}:=\operatorname{grad}_{x}\mathcal{L}(w),F_{y}:=h(x),F_{z}:=g(x)+s,F_{s}:=ZSe.$
By using these symbols and the formulation of $\nabla F(w)$ in (8), the _full_
(perturbed) Newton equation (16) defined on $T_{w}\mathcal{N}\cong
T_{x}\mathcal{M}\times\mathbb{R}^{l}\times\mathbb{R}^{m}\times\mathbb{R}^{m}$
is expanded as:

$\left(\begin{array}[]{l}\operatorname{Hess}_{x}\mathcal{L}(w)[\Delta
x]+\mathcal{H}_{x}[\Delta y]+\mathcal{G}_{x}[\Delta z]\\\
\mathcal{H}_{x}^{*}[\Delta x]\\\ \mathcal{G}_{x}^{*}[\Delta x]+\Delta s\\\
Z\Delta s+S\Delta z\end{array}\right)=\left(\begin{array}[]{l}-F_{x}\\\
-F_{y}\\\ -F_{z}\\\ -F_{s}+\mu e\end{array}\right).$ (18)

Not only does this equation contain four variables, but there is no symmetry
on the left side of the equation, so it would be unwise to solve it just like
that. By $(z,s)>0$ and the fourth line of (18), we deduce $\Delta
s=Z^{-1}\left(\mu e-F_{s}-S\Delta z\right).$ Substituting this $\Delta s$ into
the third line of (18), we get $\mathcal{G}_{x}^{*}[\Delta x]-Z^{-1}S\Delta
z=-Z^{-1}\mu e-g(x)$ and thus $\Delta
z=S^{-1}\left[Z\left(\mathcal{G}_{x}^{*}[\Delta x]+F_{z}\right)+\mu
e-F_{s}\right]$. Substituting this $\Delta z$ further into the first line of
(18) and combining it with the second line of (18) yields the following
_condensed_ Newton equation, which is defined on
$T_{x}\mathcal{M}\times\mathbb{R}^{l}$:

$\mathcal{T}(\Delta x,\Delta
y):=\left(\begin{array}[]{l}\mathcal{A}_{w}[\Delta x]+\mathcal{H}_{x}[\Delta
y]\\\ \mathcal{H}_{x}^{*}[\Delta
x]\end{array}\right)=\left(\begin{array}[]{l}c\\\ q\end{array}\right).$ (19)

where

$\displaystyle\mathcal{A}_{w}$
$\displaystyle:=\operatorname{Hess}_{x}\mathcal{L}(w)+\mathcal{G}_{x}S^{-1}Z\mathcal{G}_{x}^{*},$
(20) $\displaystyle c$
$\displaystyle:=-F_{x}-\mathcal{G}_{x}S^{-1}\left(ZF_{z}+\mu
e-F_{s}\right),q:=-F_{y}.$

Here, $c$ and $q$ are constant vectors. If we defined
$\Theta:=\mathcal{G}_{x}S^{-1}Z\mathcal{G}_{x}^{*}$, then
$\mathcal{A}_{w}=\operatorname{Hess}_{x}\mathcal{L}(w)+\Theta.$ Note that both
$\Theta$ and $\operatorname{Hess}_{x}\mathcal{L}(w)$ are operators from and to
$T_{x}\mathcal{M}.$

From the discussion above, for any $w\in\mathcal{N}$ with $(z,s)>0$, the
operator $\nabla F(w)$ in (8) is nonsingular if and only if the newly defined
operator $\mathcal{T}$ in (19) is nonsingular. Eventually, it is sufficient
for us to solve the equation (19) containing only two variables $\Delta x$ and
$\Delta y$. In fact, when we consider the case of only inequality constraints
in (RCOP), then $\Delta y$ vanishes, and only a linear equation

$\mathcal{A}_{w}[\Delta x]=c$ (21)

on $T_{x}\mathcal{M}$ needs to be solved. More importantly, the operator
$\mathcal{T}$ in the left side of (19) is symmetric, or say self-adjoint
(although often indefinite). It is trivial to check that operators $\Theta$
and $\mathcal{A}_{w}$ are self-adjoint on $T_{x}\mathcal{M}$; and thus
$\mathcal{T}$ is self-adjoint on the product vector space
$T_{x}\mathcal{M}\times\mathbb{R}^{l}$ equipped with inner product
$\left\langle(\xi_{x},\xi_{y}),(\eta_{x},\eta_{y})\right\rangle:=\left\langle\xi_{x},\eta_{x}\right\rangle_{x}+\xi_{y}^{T}\eta_{y}$.
We can also see that (19) is a saddle point problem defined on Hilbert spaces
form its special structure.

####  3.4.2 Krylov Subspace Methods on Tangent Space

Next, how to solve (19) efficiently 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 $\mathcal{A}_{w}\colon
T_{x}\mathcal{M}\to T_{x}\mathcal{M}$. Let
$d:=\operatorname{dim}T_{x}\mathcal{M}$. Unfortunately, in most cases of
practical applications, the Riemannian situation leaves us with no explicit
matrix form available for $\mathcal{A}_{w}$. This means that we can only
access $\mathcal{A}$ (subscript $w$ omitted) by inputting a vector $v$ to
return $\mathcal{A}v$.

A general approach is first to find the matrix representation
$\hat{\mathcal{A}}$ for $\mathcal{A}$ under some basis of $T_{x}\mathcal{M}$.
In detail, the full process of this approach is as follows:

  1. (Step 1)

Obtain $d$ random independent vectors on $T_{x}\mathcal{M}$.

  2. (Step 2)

Obtain an orthonormal basis $\left\\{u_{i}\right\\}_{i=1}^{d}$ of
$T_{x}\mathcal{M}$ by the modified Gram-Schmidt algorithm.

  3. (Step 3)

Compute
$[\hat{\mathcal{A}}]_{ij}:=\left\langle\mathcal{A}u_{j},u_{i}\right\rangle_{x}$
for $1\leq i\leq j\leq d$ due to symmetry, then we obtain the matrix
representation $\hat{\mathcal{A}}\in\mathbb{R}^{d\times d}$.

  4. (Step 4)

Compute $[\hat{c}]_{i}:=\left\langle c,u_{i}\right\rangle_{x}$ for $1\leq
i\leq d$, then we obtain the vector representation $\hat{c}\in\mathbb{R}^{d}$.

  5. (Step 5)

Using arbitrary linear solver to get solution $\Delta\hat{x}\in\mathbb{R}^{d}$
from matrix equation $\hat{\mathcal{A}}\Delta\hat{x}=\hat{c}$.

  6. (Step 6)

Recovery the tangent vector $\Delta x\in T_{x}\mathcal{M}$ by $\Delta
x=\sum_{i=1}^{d}(\Delta\hat{x})_{i}u_{i}$.

In Algorithm 2, at each iteration, $x$ is updated, and thus, the tangent space
$T_{x}\mathcal{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
$T_{x}\mathcal{M}$ directly. Such a method does not explicitly require a
matrix representation $\hat{\mathcal{A}}$ for $\mathcal{A}$. In general, it
only needs to call an abstract linear operator $v\mapsto\mathcal{A}v$. Since
$\mathcal{A}$ in (20) is self-adjoint but indefinite, for solving operator
equation (21), we will use the Conjugate Residual (CR) method (see [29,
ALGORITHM 6.20]) as stated in Algorithm 3.

Input: Symmetric invertible linear operator $\mathcal{A}:T_{x}\mathcal{M}\to
T_{x}\mathcal{M}$, nonzero $c\in T_{x}\mathcal{M}$ and an initial point
$v_{0}\in T_{x}\mathcal{M}$.

Output: Sequence $\left\\{v_{n}\right\\}\subset T_{x}\mathcal{M}$ such that
$\left\\{v_{n}\right\\}\to v^{*}$ and $\mathcal{A}v^{*}=c$.

Set $n\to 0$, $r_{0}:=c-Av_{0},p_{0}:=r_{0}$ and compute
$\mathcal{A}r_{0},\mathcal{A}p_{0}$ ;

while  _stopping criterion not satisfied_ do

1\. Update number $\alpha_{n}:=\langle
r_{n},\mathcal{A}r_{n}\rangle_{x}/\langle\mathcal{A}p_{n},\mathcal{A}p_{n}\rangle_{x}$
;

// Step length

2\. Set $v_{n+1}:=v_{n}+\alpha_{n}p_{n}$ ;

// Iterate point

3\. Update $r_{n+1}:=r_{n}-\alpha_{n}\mathcal{A}p_{n}$ ;

// Residual

4\. Compute $\mathcal{A}r_{n+1}$ ;

// This is the only call to $\mathcal{A}$ in while loop

5\. Update number $\beta_{n}:=\langle
r_{n+1},\mathcal{A}r_{n+1}\rangle_{x}/\langle
r_{n},\mathcal{A}r_{n}\rangle_{x}$ ;

6\. Set $p_{n+1}:=r_{n+1}+\beta_{n}p_{n}$ ;

// Conjugate direction

7\. Compute $\mathcal{A}p_{n+1}:=\mathcal{A}r_{n+1}+\beta_{n}\mathcal{A}p_{n}$
;

// No need to call $\mathcal{A}$ here

8\. $n\to n+1$;

end while

Algorithm 3 Conjugate Residual (CR) Method on Tangent Spaces for (21)

A significant feature is that the iterate points $v_{k}$, conjugate directions
$p_{k}$, and residual vectors $r_{k}:=\mathcal{A}v_{k}-c$ are all contained in
$T_{x}\mathcal{M}$. Usually, the initial point $v_{0}$ is the zero element of
$T_{x}\mathcal{M}$; the iteration terminates when the relative residual
$\left\|r_{k}\right\|/\|c\|\leq\epsilon$ for some threshold $\epsilon>0$, or
some maximum 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
$\mathcal{T}$ on $T_{x}\mathcal{M}\times\mathbb{R}^{l}$ instead of
$\mathcal{A}$ on $T_{x}\mathcal{M}$.

##  4 Preliminaries and Auxiliary Results

This section introduces the useful results that are indispensable to our
subsequent discussion.

######  Remark 4.1.

For a retraction $\mathrm{R}$ on $\mathcal{M}$ and $x\in\mathcal{M}$, by
$\mathcal{D}\mathrm{R}_{x}\left(0_{x}\right)=\mathrm{id}_{T_{x}\mathcal{M}}$
and the inverse function theorem, there exists a neighborhood $V$ of $0_{x}$
in $T_{x}\mathcal{M}$ such that $\mathrm{R}_{x}$ is a diffeomorphism on $V$;
thus, $\mathrm{R}_{x}^{-1}(y)$ is well defined for all $y\in\mathcal{M}$
sufficiently close to $x$. In this case, $\mathrm{R}_{x}(V)\subset\mathcal{M}$
is called a retractive neighborhood of $x$. Furthermore, the existence of a
_totally_ retractive neighborhood [38, Theorem 2] shows that for any
$\bar{x}\in\mathcal{M}$ there is a neighborhood $W$ of $\bar{x}$ such that
$\mathrm{R}_{x}^{-1}(y)$ is well defined for all $x,y\in W$. In what follows,
we will suppose that an appropriate neighborhood has been chosen by default
for the well-definedness of $\mathrm{R}_{x}^{-1}(y)$.

###  4.1 Vector Transport and Parallel Transport

Define the Whitney sum $T\mathcal{M}\oplus
T\mathcal{M}:=\left\\{(\eta,\xi):\eta,\xi\in
T_{x}\mathcal{M},x\in\mathcal{M}\right\\}.$ A smooth map $\mathrm{T}\colon
T\mathcal{M}\oplus T\mathcal{M}\to
T\mathcal{M}:(\eta,\xi)\mapsto\mathrm{T}_{\eta}(\xi),$ is called a _vector
transport_ on $\mathcal{M}$ if there exists an associated retraction
$\mathrm{R}$ on $\mathcal{M}$ such that $\mathrm{T}$ satisfies the following
properties for all $x\in\mathcal{M}$:

  1. 1.

Associated retraction. $\mathrm{T}_{\eta}\left(\xi\right)\in
T_{\mathrm{R}_{x}\left(\eta\right)}\mathcal{M}$ for all $\eta,\xi\in
T_{x}\mathcal{M}$.

  2. 2.

Consistency. $\mathrm{T}_{0_{x}}\left(\xi\right)=\xi$ for all $\xi\in
T_{x}\mathcal{M}$.

  3. 3.

Linearity.
$\mathrm{T}_{\eta}\left(a\xi+b\zeta\right)=a\mathrm{T}_{\eta}\left(\xi\right)+b\mathrm{T}_{\eta}\left(\zeta\right)$
for all $a,b\in\mathbb{R}$ and $\eta,\xi,\zeta\in T_{x}\mathcal{M}$.

Thus, fixing any $\eta\in T_{x}\mathcal{M}$, the map $\mathrm{T}_{\eta}\colon
T_{x}\mathcal{M}\to
T_{\mathrm{R}_{x}(\eta)}\mathcal{M}:\xi\mapsto\mathrm{T}_{\eta}(\xi),$ is a
linear operator. Additionally, $\mathrm{T}$ is _isometric_ if
$\langle\mathrm{T}_{\eta}(\xi),\mathrm{T}_{\eta}(\zeta)\rangle=\langle\xi,\zeta\rangle$
holds, for all $x\in\mathcal{M}$ and all $\xi,\zeta,\eta\in T_{x}\mathcal{M}$.
In other words, for any $\eta\in T_{x}\mathcal{M}$, the adjoint and the
inverse of $\mathrm{T}_{\eta}$ coincide, i.e.,
$\mathrm{T}_{\eta}^{*}=\mathrm{T}_{\eta}^{-1}$. There are two important
classes of vector transport as follows. Let $\mathrm{R}$ be a retraction on
$\mathcal{M}$.

  1. (1)

The _differentiated retraction_ defined by

$\mathrm{T}_{\eta}(\xi):=\mathcal{D}\mathrm{R}_{x}(\eta)[\xi],\;\eta,\xi\in
T_{x}\mathcal{M},x\in\mathcal{M},$ (22)

is a valid vector transport [1, equation (8.6)].

  2. (2)

Given a smooth curve $\gamma\colon[0,1]\to\mathcal{M}$ and
$t_{0},t_{1}\in[0,1]$, the _parallel transport_ from the tangent space at
$\gamma(t_{0})$ to the tangent space at $\gamma(t_{1})$ along $\gamma,$ is a
linear operator $\mathrm{P}_{\gamma}^{t_{1}\to t_{0}}\colon
T_{\gamma\left(t_{0}\right)}\mathcal{M}\rightarrow
T_{\gamma\left(t_{1}\right)}\mathcal{M}$ defined by
$\mathrm{P}_{\gamma}^{t_{1}\to t_{0}}(\xi)=Z\left(t_{1}\right)$, where $Z$ is
the unique parallel vector field such that $Z\left(t_{0}\right)=\xi$. Then,
for any $x\in\mathcal{M},\eta\in T_{x}\mathcal{M}$, then

$\mathrm{T}_{\eta}\left(\xi\right):=\mathrm{P}_{\gamma}^{1\to 0}(\xi)$ (23)

is a valid vector transport [1, equation (8.2)], where $\mathrm{P}_{\gamma}$
denotes the parallel transport along the curve
$t\mapsto\gamma(t):=\mathrm{R}_{x}\left(t\eta\right)$. 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 defined, e.g.,
[5, Section 10.4]. Here, we consider Lipschitz continuity with respect to a
vector transport. In what follows, let $\mathcal{M}$ be a Riemannian manifold
endowed with a vector transport $\mathrm{T}$ and an associated retraction
$\mathrm{R}$. We first consider the Lipschitz continuous gradient of scale
field $f$.

######  Definition 4.2 ([17, Definition 5.2.1]).

A function $f\colon\mathcal{M}\rightarrow\mathbb{R}$ is _Lipschitz
continuously differentiable with respect to $\mathrm{T}$ in
$\mathcal{U}\subset\mathcal{M}$_ if it is differentiable and there exists a
constant $\kappa>0$ such that, for all $x,y\in\mathcal{U}$,
$\left\|\operatorname{grad}f(y)-\mathrm{T}_{\eta}[\operatorname{grad}f(x)]\right\|\leq\kappa\|\eta\|,$
where $\eta=\mathrm{R}_{x}^{-1}y$.

Going one degree higher, let us now discuss the Lipschitz continuity of
Hessian operators. Throughout this paper, for a linear operator
$\mathcal{A}:\mathcal{E}\rightarrow\mathcal{E}^{\prime}$ between two finite-
dimensional normed vector spaces $\mathcal{E}$ and $\mathcal{E}^{\prime}$, the
(operator) norm of $\mathcal{A}$ is defined by
$\|\mathcal{A}\|:=\sup\left\\{\|\mathcal{A}v\|_{\mathcal{E}^{\prime}}:v\in\mathcal{E},\|v\|_{\mathcal{E}}=1\right.$,
or, $\left.\|v\|_{\mathcal{E}}\leq 1\right\\}.$ The inverse of
$\mathrm{T}_{\eta}$ is needed in the following definitions, so we can assume
that vector transport $\mathrm{T}$ is isometric, e.g., parallel transport in
(23). In fact, there are many ways to construct isometric vector transports;
see [18, Section 2.3].

######  Definition 4.3 ([18, Assumption 3]).

A function $f\colon\mathcal{M}\rightarrow\mathbb{R}$ is _twice Lipschitz
continuously differentiable with respect to $\mathrm{T}$ in
$\mathcal{U}\subset\mathcal{M}$_ if it is twice differentiable and there
exists a constant $\kappa>0$ such that, for all $x,y\in\mathcal{U}$,
$\|\operatorname{Hess}f(y)-\mathrm{T}_{\eta}\operatorname{Hess}f(x)\mathrm{T}_{\eta}^{-1}\|\leq\kappa
d(x,y),$ where $\eta=\mathrm{R}_{x}^{-1}y$.

######  Lemma 4.4 ([18, Lemma 4]).

If $f\colon\mathcal{M}\rightarrow\mathbb{R}$ is $C^{3}$, then for any
$\bar{x}\in\mathcal{M}$ and any isometric vector transport $\mathrm{T}$, there
exists a neighborhood $\mathcal{U}$ of $\bar{x}$ such that $f$ is twice
Lipschitz continuously differentiable with respect to $\mathrm{T}$ in
$\mathcal{U}$.

If the operator, $\operatorname{Hess}f(x)$, above is replaced by a general
covariant derivative $\nabla F(x)$, we can get the next results in a similar
way. Lemma 4.6 can be proven in the same way as Lemma 4.4.

######  Definition 4.5.

Given a vector field $F$ on $\mathcal{M}$. The map $x\mapsto\nabla F(x)$ is
_Lipschitz continuous with respect to $\mathrm{T}$ in
$\mathcal{U}\subset\mathcal{M}$_ if there exists a constant $\kappa>0$ such
that, for all $x,y\in\mathcal{U}$, it holds that $\|\nabla
F(y)-\mathrm{T}_{\eta}\nabla F(x)\mathrm{T}_{\eta}^{-1}\|\leq\kappa d(x,y),$
where $\eta=\mathrm{R}_{x}^{-1}y$.

######  Lemma 4.6.

If $F$ is a $C^{2}$ vector field, then for any $\bar{x}\in\mathcal{M}$ and any
isometric vector transport $\mathrm{T}$, there exists a neighborhood
$\mathcal{U}$ of $\bar{x}$ such that the map $x\mapsto\nabla F(x)$ is
Lipschitz continuous with respect to $\mathrm{T}$ in $\mathcal{U}$.

###  4.3 Auxiliary Lemmas

Notice that in the previous subsection on the definitions of Lipschitz
continuity, we used $\|\eta\|$ with $\eta=\mathrm{R}_{x}^{-1}y$ or $d(x,y)$ to
denote the upper bound on the right-hand side. The next lemma shows that the
two are not essentially different. When $\mathcal{M}=\mathbb{R}^{n}$, both
reduce to $\|x-y\|$.

######  Lemma 4.7 ([18, Lemma 2]).

Let $\mathcal{M}$ be a Riemannian manifold with a retraction $\mathrm{R}$ and
let $\bar{x}\in\mathcal{M}$. Then,

  1. (i)

there exist $a_{0},a_{1},\delta_{a_{0},a_{1}}>0$ such that for all $x$ in a
sufficiently small neighborhood of $\bar{x}$ and all $\xi,\eta\in
T_{x}\mathcal{M}$ with $\|\xi\|,\|\eta\|\leq\delta_{a_{0},a_{1}}$, one has
$a_{0}\|\xi-\eta\|\leq d(\mathrm{R}_{x}(\eta),\mathrm{R}_{x}(\xi))\leq
a_{1}\|\xi-\eta\|.$ In particular, $a_{0}\|\xi\|\leq
d(x,\mathrm{R}_{x}(\xi))\leq a_{1}\|\xi\|$ when $\eta=0$;

  2. (ii)

there exist $a_{0},a_{1}>0$ such that for all $x$ in a sufficiently small
neighborhood of $\bar{x}$, one has $a_{0}\|\xi\|\leq d(x,\bar{x})\leq
a_{1}\|\xi\|\text{ where }\xi=\mathrm{R}_{\bar{x}}^{-1}(x).$

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 field and $\bar{x}\in\mathcal{M}$. Then there
exist a neighborhood $\mathcal{U}$ of $\bar{x}$ and a constant $c_{1}\geq 0$
such that for all $x,y\in\mathcal{U}$,

$\left\|\mathrm{P}_{\gamma}^{0\to
1}[F(y)]-F(x)-\int_{0}^{1}\mathrm{P}_{\gamma}^{0\to t}\nabla
F(\gamma(t))\mathrm{P}_{\gamma}^{t\to 0}[\eta]\mathrm{d}t\right\|\leq
c_{1}\|\eta\|^{2},$

where $\eta=\mathrm{R}_{x}^{-1}(y)$ and $\mathrm{P}_{\gamma}$ is the parallel
transport along the curve $\gamma(t):=\mathrm{R}_{x}\left(t\eta\right)$.
Moreover, if $\mathrm{R}=\operatorname{Exp}$, then indeed $c_{1}=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 field and $\bar{x}\in\mathcal{M}$. Then there
exist a neighborhood $\mathcal{U}$ of $\bar{x}$ and a constant $c_{2}>0$ such
that for all $x\in\mathcal{U}$,

$\left\|\mathrm{P}_{\gamma}^{0\to 1}[F(x)]-F(\bar{x})-\nabla
F(\bar{x})[\eta]\right\|\leq c_{2}d^{2}(\bar{x},x),$

where $\eta=\mathrm{R}_{\bar{x}}^{-1}x$ and $\mathrm{P}_{\gamma}$ is the
parallel transport along the curve $\gamma(t):=\mathrm{R}_{\bar{x}}(t\eta)$.

###### Proof.

Let $LHS:=\|\mathrm{P}_{\gamma}^{0\to 1}[F(x)]-F(\bar{x})-\nabla
F(\bar{x})[\eta]\|.$ It follows that

$\displaystyle LHS$ $\displaystyle\leq\left\|\mathrm{P}_{\gamma}^{0\to
1}[F(x)]-F(\bar{x})-\int_{0}^{1}\mathrm{P}_{\gamma}^{0\to t}\nabla
F(\gamma(t))\mathrm{P}_{\gamma}^{t\to 0}[\eta]\mathrm{d}t\right\|$
$\displaystyle\quad+\left\|\int_{0}^{1}\mathrm{P}_{\gamma}^{0\to t}\nabla
F(\gamma(t))\mathrm{P}_{\gamma}^{t\to 0}[\eta]\mathrm{d}t-\nabla
F(\bar{x})[\eta]\right\|$ $\displaystyle\leq
c_{1}\|\eta\|^{2}+\Big{\|}\int_{0}^{1}\left(\mathrm{P}_{\gamma}^{0\to t}\nabla
F(\gamma(t))\mathrm{P}_{\gamma}^{t\to 0}-\nabla
F(\bar{x})\right)[\eta]\mathrm{d}t\Big{\|}.\text{ (by Lemma
\ref{lem:FundamentalTheoremOfCalculusRiemannian})}$

Let $\theta:=\left\|\int_{0}^{1}\left(\mathrm{P}_{\gamma}^{0\to t}\nabla
F(\gamma(t))\mathrm{P}_{\gamma}^{t\to 0}-\nabla
F(\bar{x})\right)[\eta]\mathrm{d}t\right\|$. Note that

$\displaystyle\theta$
$\displaystyle\leq\int_{0}^{1}\left\|\mathrm{P}_{\gamma}^{0\to t}\nabla
F(\gamma(t))\mathrm{P}_{\gamma}^{t\to 0}-\nabla
F(\bar{x})\right\|\left\|\eta\right\|\mathrm{d}t$
$\displaystyle\leq\int_{0}^{1}c_{0}d(\bar{x},\mathrm{R}_{\bar{x}}(t\eta))\left\|\eta\right\|\mathrm{d}t\text{
(by Lemma \ref{lem:LipCountNabulaF})}$
$\displaystyle\leq\int_{0}^{1}c_{0}a_{1}t\|\eta\|\left\|\eta\right\|\mathrm{d}t=\frac{1}{2}c_{0}a_{1}\|\eta\|^{2}.\text{
(by (i) of Lemma \ref{lem:distR1})}$

Combining the above results yields

$LHS\leq(c_{1}+\frac{1}{2}c_{0}a_{1})\|\eta\|^{2}\leq(c_{1}+\frac{1}{2}c_{0}a_{1})/a_{0}^{2}d^{2}(\bar{x},x),$

where the last inequality comes from (ii) of Lemma 4.7. Letting
$c_{2}:=(c_{1}+\frac{1}{2}c_{0}a_{1})/a_{0}^{2}$ completes the proof. ∎

We end this section with the following useful lemmas.

######  Lemma 4.10 ([12, Lemma 3.2]).

Given a vector field $F$ on $\mathcal{M}$. If the map $p\mapsto\nabla F(p)$ is
continuous at $p^{*}$ and $\nabla F(p^{*})$ is nonsingular, then there exist a
neighborhood $\mathcal{U}$ of $p^{*}$ and a constant $\Xi>0$ such that, for
all $p\in\mathcal{U}$, $\nabla F(p)$ is nonsingular and $\left\|\nabla
F(p)^{-1}\right\|\leq\Xi$.

######  Lemma 4.11 ([14, Lemma 14.5]).

Let $F$ be a $C^{2}$ vector field on $\mathcal{M}$ and $p^{*}\in\mathcal{M}$.
If $F(p^{*})=0$ and $\nabla F(p^{*})$ is nonsingular, then there exist a
neighborhood $\mathcal{U}$ of $p^{*}$ and constants $c_{3},c_{4}>0$ such that,
for all $p\in\mathcal{U}$, $c_{3}d(p,p^{*})\leq\|F(p)\|\leq c_{4}d(p,p^{*}).$

######  Lemma 4.12 ([7, Lemma 3.5]).

Let $u\in T_{p}\mathcal{M}$ such that $\operatorname{Exp}_{p}(u)$ exists and
$v\in T_{p}\mathcal{M}\cong T_{p}\left(T_{p}\mathcal{M}\right)$. Then
$\langle\mathcal{D}\operatorname{Exp}_{p}(u)[u],\mathcal{D}\operatorname{Exp}_{p}(u)[v]\rangle=\langle
u,v\rangle.$ In particular, $\left\|\mathcal{D}\operatorname{Exp}_{p}(\lambda
u)[u]\right\|=\left\|u\right\|$ holds all $\lambda\geq 0.$

##  5 Local Convergence

Here, for any two nonnegative sequences $\left\\{u_{k}\right\\}$ and
$\left\\{v_{k}\right\\}$, we write $u_{k}=O(v_{k})$ if there is a constant
$M>0$ such that $u_{k}\leq Mv_{k}$ for all sufficiently large $k$; and we
write $u_{k}=o(v_{k})$ if $v_{k}>0$ and the sequence of ratios
$\\{u_{k}/v_{k}\\}$ approaches zero. In this section, we will establish local
convergence of our 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 for solving the singularity problem (2), which can be stated as
Algorithm 4.

Input: A vector field $F\in\mathfrak{X}(\mathcal{M})$, an initial point
$p_{0}\in\mathcal{M}$ and a retraction $\mathrm{R}$ on $\mathcal{M}$. Set
$\mu_{0}>0$.

Output: Sequence $\\{p_{k}\\}\subset\mathcal{M}$ such that
$\left\\{p_{k}\right\\}\rightarrow p^{*}$ and $F(p^{*})=0_{p^{*}}$.

Set $k\to 0$;

while  _Stopping criterion not satisfied_ do

1\. Obtain $\xi_{k}\in T_{p_{k}}\mathcal{M}$ by solving the perturbed Newton
equation:  $\nabla F(p_{k})[\xi_{k}]=-F(p_{k})+\mu_{k}\hat{e};$ (24)

2\. Choose a (damped) step size $0<\alpha_{k}\leq 1$;

3\. Compute the next point as
$p_{k+1}:=\mathrm{R}_{p_{k}}(\alpha_{k}\xi_{k})$;

4\. Choose $0<\mu_{k+1}<\mu_{k}$;

5\. $k\to k+1$;

end while

Algorithm 4 Perturbed Damped Riemannian Newton Method for (2)

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 $\mu_{k}\hat{e}$, while “damped” means using $\alpha_{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:

  1. (B1)

There exists $p^{*}\in\mathcal{M}$ such that $F(p^{*})=0_{p^{*}}$;

  2. (B2)

The covariant derivative $\nabla F(p^{*})$ is nonsingular;

  3. (B3)

The vector field $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 $\mu_{k}$ and $\alpha_{k}$
according to the two schemes that Proposition 5.1 gives. We can see that
either scheme will have $\mu_{k}\to 0$ and $\alpha_{k}\to 1$, which makes
Algorithm 4 eventually reduce to Algorithm 1 when $k$ is sufficiently 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
$\mu_{k},\alpha_{k}$ as follows; then there exists a constant $\delta>0$ such
that for all $p_{0}\in\mathcal{M}$ with $d(p_{0},p^{*})<\delta,$ the sequence
$\left\\{p_{k}\right\\}$ is well defined. Furthermore,

  1. (1)

if we choose $\mu_{k}=o(\|F(p_{k})\|)$ and $\alpha_{k}\to 1$, then $p_{k}\to
p^{*}$ superlinearly;

  2. (2)

if we choose $\mu_{k}=O(\|F(p_{k})\|^{2})$ and $1-\alpha_{k}=O(\|F(p_{k})\|)$,
then $p_{k}\to p^{*}$ quadratically.

###### Proof.

By (B2)-(B3), Lemma 4.6 and Lemma 4.10, we can let $p_{k}$ be sufficiently
close to $p^{*}$ such that $\nabla F(p_{k})$ is nonsingular, and
$\left\|\nabla F(p_{k})^{-1}\right\|\leq\Xi$ for some constant $\Xi$. Then,
the next iterate point,

$p_{k+1}:=\mathrm{R}_{p_{k}}[\alpha_{k}\nabla
F(p_{k})^{-1}(-F(p_{k})+\mu_{k}\hat{e})],$

is well defined in Algorithm 4, and it follows from
$p^{*}=\mathrm{R}_{p_{k}}(\eta)$ with $\eta:=\mathrm{R}_{p_{k}}^{-1}p^{*}$ and
(i) of Lemma 4.7 that

$\displaystyle d(p_{k+1},p^{*})$ $\displaystyle\leq
a_{1}\|\eta-\alpha_{k}\nabla F(p_{k})^{-1}(-F(p_{k})+\mu_{k}\hat{e})\|$
$\displaystyle=a_{1}\|\eta+\alpha_{k}\nabla
F(p_{k})^{-1}(F(p_{k})-\mu_{k}\hat{e})\|.$ (25)

Let $r_{k}:=\eta+\alpha_{k}\nabla F(p_{k})^{-1}(F(p_{k})-\mu_{k}\hat{e})$.
Algebraic manipulations show that

$r_{k}=(1-\alpha_{k})\eta+\alpha_{k}\nabla F(p_{k})^{-1}[\nabla
F(p_{k})\eta+F(p_{k})-\mathrm{P}_{\gamma}^{0\to 1}F(p^{*})-\mu_{k}\hat{e}],$

where $\mathrm{P}_{\gamma}$ is the parallel transport along the curve
$\gamma(t)=\mathrm{R}_{p_{k}}(t\eta)$ and $F(p^{*})=0.$ Thus, using
$\left\|\eta\right\|\leq d(p_{k},p^{*})/a_{0}$ from (ii) of Lemma 4.7 and
Lemma 4.9, we have

$\displaystyle\|r_{k}\|$
$\displaystyle\leq(1-\alpha_{k})\|\eta\|+\alpha_{k}\|\nabla
F(p_{k})^{-1}\|\|\mathrm{P}_{\gamma}^{0\to 1}F(p^{*})-F(p_{k})-\nabla
F(p_{k})\eta\|$ $\displaystyle\quad+\alpha_{k}\|\nabla
F(p_{k})^{-1}\|\|\hat{e}\|\mu_{k}$
$\displaystyle\leq(1-\alpha_{k})d(p_{k},p^{*})/a_{0}+\alpha_{k}\|\nabla
F(p_{k})^{-1}\|c_{2}d^{2}(p_{k},p^{*})$ $\displaystyle\quad+\alpha_{k}\|\nabla
F(p_{k})^{-1}\|\|\hat{e}\|\mu_{k}$
$\displaystyle\leq(1-\alpha_{k})d(p_{k},p^{*})/a_{0}+\Xi
c_{2}d^{2}(p_{k},p^{*})+\Xi\|\hat{e}\|\mu_{k}.$

Combining the above with (25), we conclude that

$d(p_{k+1},p^{*})\leq\kappa_{1}(1-\alpha_{k})d(p_{k},p^{*})+\kappa_{2}d^{2}(p_{k},p^{*})+\kappa_{3}\mu_{k}$
(26)

for some positive constants $\kappa_{1},\kappa_{2},\kappa_{3}$. On the other
hand, by Lemma 4.11, we have

$\left\|F(p_{k})\right\|=O(d(p_{k},p^{*})).$ (27)

In what follows, we prove assertions (1) and (2).

(1) Suppose that $\alpha_{k}\to 1$ and $\mu_{k}=o(\|F(p_{k})\|)$, which
together with (27) imply $\mu_{k}=o(d(p_{k},p^{*}))$. By (26), we have

$\frac{d(p_{k+1},p^{*})}{d(p_{k},p^{*})}\leq\kappa_{1}(1-\alpha_{k})+\kappa_{2}d(p_{k},p^{*})+\kappa_{3}\frac{\mu_{k}}{d(p_{k},p^{*})}.$
(28)

We can take $\delta$ sufficiently small and $k$ sufficiently large, if
necessary, to conclude that
$d(p_{k+1},p^{*})<\frac{1}{2}d(p_{k},p^{*})<\delta.$ Thus, $p_{k+1}\in
B_{\delta}(p^{*})$, the open ball of radius $\delta$ centered at $p^{*}$ on
$\mathcal{M}$. By induction, it is easy to show that the sequence
$\left\\{p_{k}\right\\}$ is well defined and converges to $p^{*}$. Taking the
limit of both sides of (28) proves superlinear convergence.

(2) Again, we start from (26) and rewrite it as:

$\displaystyle
d(p_{k+1},p^{*})=(1-\alpha_{k})O(d(p_{k},p^{*}))+O(d^{2}(p_{k},p^{*}))+O(\mu_{k}).$
(29)

Suppose that $1-\alpha_{k}=O(\left\|F(p_{k})\right\|)$ and
$\mu_{k}=O(\left\|F(p_{k})\right\|^{2})$. Using (27), the above reduces to
$d(p_{k+1},p^{*})=O(d^{2}(p_{k},p^{*})).$ This implies that there exists a
constant $\nu$ such that $d(p_{k+1},p^{*})\leq\nu d^{2}(p_{k},p^{*})$, and
hence, $d(p_{k+1},p^{*})\leq\nu d^{2}(p_{k},p^{*})\leq\nu\delta^{2}<\delta,$
if $\delta$ is sufficiently small. Again, by induction,
$\left\\{p_{k}\right\\}$ converges to $p^{*}$ quadratically. ∎

###  5.2 Local Convergence of Algorithm 2

Next, lemma shows the relationship between the parameter $\gamma_{k}$ and step
size $\alpha_{k}$ in Algorithm 2.

######  Lemma 5.2.

Consider the Algorithm 2 for solving the problem (RCOP). Let (A1) and (A3)
hold at some $w^{*}=(x^{*},y^{*},z^{*},s^{*})$ and $\alpha_{k}$ be as in (17).
Define a constant,

$\Pi:=2\max\left\\{\max_{i}\left\\{1/(s^{*})_{i}\mid(s^{*})_{i}>0\right\\},\max_{i}\left\\{1/(z^{*})_{i}\mid(z^{*})_{i}>0\right\\}\right\\}.$

For $\gamma_{k}\in(0,1)$, if $\Pi\left\|\Delta w_{k}\right\|\leq\gamma_{k},$
then $0\leq 1-\alpha_{k}\leq(1-\gamma_{k})+\Pi\left\|\Delta w_{k}\right\|.$

###### Proof.

Notice that the fourth line of (18) yields

$S_{k}^{-1}\Delta s_{k}+Z_{k}^{-1}\Delta z_{k}=\mu_{k}(S_{k}Z_{k})^{-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 Riemannian case. ∎

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
$\gamma_{k}$.

######  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 $\mu_{k},\gamma_{k}$ as follows; then there
exists a constant $\delta>0$ such that, for all $w_{0}\in\mathcal{N}$ with
$d(w_{0},w^{*})<\delta,$ the sequence $\left\\{w_{k}\right\\}$ is well
defined. Furthermore,

  1. (1)

if we choose $\mu_{k}=o(\left\|F(w_{k})\right\|)$ and $\gamma_{k}\to 1$, then
$w_{k}\to w^{*}$ superlinearly;

  2. (2)

if we choose $\mu_{k}=O(\left\|F(w_{k})\right\|^{2})$ and
$1-\gamma_{k}=O(\left\|F(w_{k})\right\|)$, then $w_{k}\to w^{*}$
quadratically.

###### Proof.

We only prove (2) because (1) can be proven in the same way. Let $w_{k}$ be
such that $d(w_{k},w^{*})<\delta$ for sufficiently small $\delta$. From
Proposition 3.1, (A1)-(A4) shows that the KKT vector field $F$ satisfies
(B1)-(B3); thus, the discussion in the proof of Proposition 5.1 applies to KKT
vector field $F$ as well, simply by replacing the symbol $p$ with $w$. Since
we choose $\mu_{k}=O(\|F(w_{k})\|^{2})$, and $\|F(w_{k})\|=O(d(w_{k},w^{*}))$
by (27), we obtain $\mu_{k}=O(d^{2}(w_{k},w^{*})).$ Thus,

$\displaystyle\left\|\Delta w_{k}\right\|$ $\displaystyle=\left\|\nabla
F(w_{k})^{-1}(-F(w_{k})+\mu_{k}\hat{e})\right\|\text{ (by
(\ref{eq:prototype}))}$
$\displaystyle\leq\Xi(\left\|F(w_{k})\right\|+\mu_{k}\|\hat{e}\|)\text{ (by
Lemma \ref{lem:nonsingular})}$ $\displaystyle\leq O(\|F(w_{k})\|)+O(\mu_{k})$
$\displaystyle=O(d(w_{k},w^{*}))+O(d^{2}(w_{k},w^{*}))=O(d(w_{k},w^{*})).$

Since $\delta$ is sufficiently small, from the above inequalities, the
condition of Lemma 5.2 is satisfied, i.e., $\left\|\Delta
w_{k}\right\|\leq\hat{\gamma}/\Pi\leq\gamma_{k}/\Pi$ for a constant
$\hat{\gamma}\in(0,1).$ Hence, by
$1-\gamma_{k}=O(\left\|F\left(w_{k}\right)\right\|)$, one has
$1-\alpha_{k}\leq(1-\gamma_{k})+\Pi\left\|\Delta
w_{k}\right\|=(1-\gamma_{k})+O(d(w_{k},w^{*}))=O(d(w_{k},w^{*})).$ Finally,
from (29), we have
$d(w_{k+1},w^{*})=(1-\alpha_{k})O(d(w_{k},w^{*}))+O(d^{2}(w_{k},w^{*}))+O(\mu_{k})=O(d^{2}(w_{k},w^{*})).$
This completes the proof. ∎

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 online 333See
…/LocalRIPM/PrototypeRIPM.m in https://doi.org/10.5281/zenodo.10612799, where
from (2) of Theorem 5.3, in practice we set parameters
$\mu_{k+1}=\min\\{\mu_{k}/1.5,0.5\left\|F(w_{k})\right\|^{2}\\};$ and
$\gamma_{k+1}=\max\\{\hat{\gamma},1-\left\|F(w_{k})\right\|\\}$ with
$\hat{\gamma}=0.5$.

##  6 Global Algorithm

The globally convergent version of our RIPM uses the classical line search
described in [11]. The following considerations and definitions are needed in
order to describe it compactly. For simplicity, we often omit the subscript of
iteration count $k$.

Given the current point $w=(x,y,z,s)$ and $\Delta w=(\Delta x,\Delta y,\Delta
z,\Delta s)$, the next iterate is obtained along a curve on product manifold
$\mathcal{N}$, i.e., $\alpha\mapsto
w(\alpha):=\bar{\mathrm{R}}_{w}(\alpha\Delta w)$ with some step size
$\alpha>0$, see (6) for $\bar{\mathrm{R}}_{w}$. By introducing
$w(\alpha)=(x(\alpha),y(\alpha),z(\alpha),s(\alpha))$, we have
$x(\alpha)=\mathrm{R}_{x}(\alpha\Delta x)$, $y(\alpha)=y+\alpha\Delta y$,
$z(\alpha)=z+\alpha\Delta z$, and $s(\alpha)=s+\alpha\Delta s$. For a given
starting point $w_{0}\in\mathcal{N}$ with $(z_{0},s_{0})>0$, let us set two
constants

$\tau_{1}:=\min(Z_{0}S_{0}e)/(z_{0}^{T}s_{0}/m),\;\tau_{2}:=z_{0}^{T}s_{0}/\left\|F(w_{0})\right\|.$

As well, define two functions
$f^{I}(\alpha):=\min(Z(\alpha)S(\alpha)e)-\gamma\tau_{1}z(\alpha)^{T}s(\alpha)/m$,
$f^{II}(\alpha):=z(\alpha)^{T}s(\alpha)-\gamma\tau_{2}\|F(w(\alpha))\|,$ where
$\gamma\in(0,1)$ is a constant. For $i=I,II$, define

$\alpha^{i}:=\max_{\alpha\in(0,1]}\left\\{\alpha:f^{i}(t)\geq 0,\text{ for all
}t\in(0,\alpha]\right\\},$ (30)

i.e., $\alpha^{i}$ are either one or the smallest positive root for the
functions $f^{i}(\alpha)$ in $(0,1]$.

Moreover, we define the merit function $\varphi\colon\mathcal{N}\to\mathbb{R}$
by $\varphi(w):=\|F(w)\|^{2}$; accordingly, we have
$\operatorname{grad}\varphi(w)=2\nabla F(w)^{*}[F(w)],$ where symbol $*$ means
its adjoint operator. Let $\|\cdot\|_{1}$, $\|\cdot\|_{2}$ be $l_{1},l_{2}$
vector norms. Note that
$\|F(w)\|_{w}^{2}=\left\|\operatorname{grad}_{x}\mathcal{L}(w)\right\|_{x}^{2}+\|h(x)\|_{2}^{2}+\|g(x)+s\|_{2}^{2}+\|ZSe\|_{2}^{2}$
by (5). Moreover, for any nonnegative $z,s\in\mathbb{R}^{m}$, one has
$\left\|ZSe\right\|_{2}\leq
z^{T}s=\|ZSe\|_{1}\leq\sqrt{m}\left\|ZSe\right\|_{2}.$ Hence,

$\left\|ZSe\right\|_{2}/\sqrt{m}\leq
z^{T}s/\sqrt{m}\leq\left\|ZSe\right\|_{2}\leq\|F(w)\|.$ (31)

Now, using the above definitions, the globally convergent RIPM can be stated
as Algorithm 5.

Input: An initial point
$w_{0}=\left(x_{0},y_{0},z_{0},s_{0}\right)\in\mathcal{N}$ with
$\left(z_{0},s_{0}\right)>0$ and a retraction $\mathrm{R}$ on $\mathcal{M}$.
$\theta\in(0,1),\beta\in(0,0.5],\gamma_{-1}\in(0.5,1).$

Output: Sequence $\left\\{w_{k}\right\\}\subset\mathcal{N}$ such that
$\left\\{w_{k}\right\\}\rightarrow w^{*}$ and $w^{*}$ satisfies the KKT
conditions (4).

Set $k\to 0$;

while  _Stopping criterion not satisfied_ do

1\. Set
$\sigma_{k}\in(0,1),\rho_{k}\in[z_{k}^{T}s_{k}/m,\left\|F(w_{k})\right\|/\sqrt{m}]$;

2\. Obtain $\Delta w_{k}=\left(\Delta x_{k},\Delta y_{k},\Delta z_{k},\Delta
s_{k}\right)\in T_{w_{k}}\mathcal{N}$ by solving the following linear
equation:  $\nabla F(w_{k})[\Delta
w_{k}]=-F(w_{k})+\sigma_{k}\rho_{k}\hat{e};$ (32)

3\. Step size selection:

(3a) Centrality condition: Set $\gamma_{k}\in(0.5,\gamma_{k-1})$ and
$\bar{\alpha}_{k}=\min\\{\alpha_{k}^{I},\alpha_{k}^{II}\\}$ from (30);

(3b) Sufficient decrease condition: Let $\alpha_{k}:=\bar{\alpha}_{k}$;

while  _$\alpha_{k}$ does not satisfy the condition
$\varphi(\bar{\mathrm{R}}_{w_{k}}(\alpha_{k}\Delta
w_{k}))-\varphi(w_{k})\leq\alpha_{k}\beta\left\langle\operatorname{grad}\varphi_{k},\Delta
w_{k}\right\rangle$_ do

$\alpha_{k}:=\theta\alpha_{k}$;

end while

4\. Compute the next point as
$w_{k+1}:=\bar{\mathrm{R}}_{w_{k}}(\alpha_{k}\Delta w_{k})$;

5\. $k\to k+1$;

end while

Algorithm 5 Global Convergent Algorithm of RIPM for (RCOP)

Compared to the prototype Algorithm 2, the global Algorithm 5 involves a more
elaborate choice of step size. Regarding the centrality condition (3a), it
does not differ in any way from the Euclidean setting; thus, references [11,
3, 9] show that $\bar{\alpha}_{k}$ is well-defined thereby ensuring that
$z_{k}$, $s_{k}$ are positive. In the following, we focus on the sufficient
decrease (Armijo) condition (3b):

$\varphi(\bar{\mathrm{R}}_{w_{k}}(\alpha_{k}\Delta
w_{k}))-\varphi(w_{k})\leq\alpha_{k}\beta\left\langle\operatorname{grad}\varphi_{k},\Delta
w_{k}\right\rangle,$ (33)

where
$\operatorname{grad}\varphi_{k}\equiv\operatorname{grad}\varphi\left(w_{k}\right)$
for short. With a slight abuse of notation $\varphi$, at the current point $w$
and direction $\Delta w$, we define a real-to-real function
$\alpha\mapsto\varphi(\alpha):=\varphi\left(\bar{\mathrm{R}}_{w}(\alpha\Delta
w)\right)$; then, it follows from the definition of a retraction and the chain
rule that

$\varphi^{\prime}(0)=\mathcal{D}\varphi(\bar{\mathrm{R}}_{w}(0))\left[\mathcal{D}\bar{\mathrm{R}}_{w}(0)[\Delta
w]\right]=\mathcal{D}\varphi(w)[\Delta
w]=\langle\operatorname{grad}\varphi(w),\Delta w\rangle.$

Hence, $\varphi_{k}^{\prime}(0)=\langle\operatorname{grad}\varphi_{k},\Delta
w_{k}\rangle$ at the $k$-th iterate $w_{k}$; and then the Armijo condition
(33) reads
$\varphi_{k}(\alpha_{k})-\varphi_{k}(0)\leq\alpha_{k}\beta\varphi_{k}^{\prime}(0)$
as usual. If $\varphi_{k}^{\prime}(0)<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 $\Delta w_{k}$ generated by (32)
ensures the descent of the merit function.

######  Lemma 6.1.

If the direction $\Delta w_{k}$ is the solution of equation (32), then

$\langle\operatorname{grad}\varphi(w_{k}),\Delta
w_{k}\rangle=2(-\|F(w_{k})\|^{2}+\sigma_{k}\rho_{k}z_{k}^{T}s_{k}).$

In this case, $\Delta w_{k}$ is a descent direction for $\varphi(w)$ at
$w_{k}$ if and only if $\rho_{k}<\|F(w_{k})\|^{2}/\sigma_{k}z_{k}^{T}s_{k}.$

###### Proof.

The iteration count $k$ is omitted. Let $\Delta w$ be given by (32). Then, we
have $\langle\operatorname{grad}\varphi(w),\Delta w\rangle=\left\langle
2\nabla F(w)^{*}[F(w)],\Delta w\right\rangle=2\langle F(w),\nabla F(w)[\Delta
w]\rangle=2\langle F(w),-F(w)+\sigma\rho\hat{e}\rangle=2(-\langle
F(w),F(w)\rangle+\sigma\rho\langle F(w),\hat{e}\rangle).$ Then, by definition
of $\hat{e}$ in (13), $\langle F(w),\hat{e}\rangle=\langle
ZSe,e\rangle=z^{T}s$ completes the proof. ∎

The next proposition shows that Algorithm 5 can generate the monotonically
nonincreasing sequence $\left\\{\varphi_{k}\right\\}$. Note that
$\varphi_{k+1}\equiv\varphi_{k}\left(\alpha_{k}\right)$ and
$\varphi_{k}\equiv\varphi_{k}(0)\equiv\varphi\left(w_{k}\right)$.

######  Proposition 6.2.

If $\|F(w_{k})\|\neq 0$, then the direction $\Delta w_{k}$ generated by
Algorithm 5 is a descent direction for $\varphi(w)$ at $w_{k}$. Moreover, if
the Armijo condition (33) is satisfied, then

$\varphi_{k}(\alpha_{k})\leq\left[1-2\alpha_{k}\beta(1-\sigma_{k})\right]\varphi_{k}(0).$

Thus, the sequence $\left\\{\varphi_{k}\right\\}$ is monotonically
nonincreasing.

###### Proof.

The iteration count $k$ is omitted. Suppose that
$\rho\leq\left\|F(w)\right\|/\sqrt{m}$ and $\Delta w$ is given by (32), we
have

$\displaystyle\varphi^{\prime}(0)$
$\displaystyle=\langle\operatorname{grad}\varphi(w),\Delta
w\rangle=2(-\varphi(w)+\sigma\rho z^{T}s)\text{ (by Lemma \ref{lem:dr})}$
$\displaystyle\leq 2(-\varphi(w)+\sigma\left\|F(w)\right\|z^{T}s/\sqrt{m})$
$\displaystyle\leq 2(-\varphi(w)+\sigma\left\|F(w)\right\|^{2})\text{ (by
(\ref{genlma2}))}$ $\displaystyle=-2(1-\sigma)\varphi(w)<0.$ (34)

Thus, in Algorithm 5, $\Delta w$ is a descent direction for the merit function
$\varphi$ at $w$. Alternatively, by Lemma 6.1, it is sufficient to show that
$\|F(w)\|/\sqrt{m}<\|F(w)\|^{2}/\sigma z^{T}s.$ By $\sigma
z^{T}s<z^{T}s\leq\sqrt{m}\|F(w)\|;$ then, $1/\sqrt{m}<\|F(w)\|/\sigma z^{T}s.$
Multiplying both sides by $\|F(w)\|$ yields the result.

Moreover, if condition (33) is satisfied, then by (34), we have
$\varphi(\alpha)\leq\varphi(0)+\alpha\beta\langle\operatorname{grad}\varphi(w),\Delta
w\rangle\leq\varphi(0)+\alpha\beta(-2(1-\sigma)\varphi(0))=[1-2\alpha\beta(1-\sigma)]\varphi(0).$
Note that in Algorithm 5, we set $\beta\in(0,1/2],\sigma\in(0,1)$, and
$\alpha\in(0,1]$, which imply that the sequence $\left\\{\varphi_{k}\right\\}$
is monotonically nonincreasing. ∎

Finally, we need to make the following assumptions. For $\epsilon\geq 0$, let

$\Omega(\epsilon):=\left\\{w\in\mathcal{N}\mid\epsilon\leq\varphi(w)\leq\varphi_{0},\min(ZSe)/(z^{T}s/m)\geq\tau_{1}/2,z^{T}s/\|F(w)\|\geq\tau_{2}/2\right\\}.$

  1. (C1)

In the set $\Omega(0)$, $f,h,$ and $g$ are smooth functions;
$\left\\{\operatorname{grad}h_{i}(x)\right\\}_{i=1}^{l}$ is linearly
independent for all $x$; and the map $w\mapsto\nabla F(w)$ is Lipschitz
continuous (with respect to parallel transport);

  2. (C2)

The sequences $\left\\{x_{k}\right\\}$ and $\left\\{z_{k}\right\\}$ are
bounded [10, 3];

  3. (C3)

In any compact subset of $\Omega(0)$, $\nabla F(w)$ is nonsingular.

Given the above assumptions, we can now prove the following statement.

######  Theorem 6.3 (Global Convergence of RIPM).

Let $\left\\{w_{k}\right\\}$ be generated by Algorithm 5 with
$\mathrm{R}=\operatorname{Exp}$ and $\left\\{\sigma_{k}\right\\}\subset(0,1)$
be bounded away from zero and one. Let $\varphi$ be Lipschitz continuous on
$\Omega(0)$. If (C1)-(C3) hold, then
$\left\\{\left\|F(w_{k})\right\|\right\\}$ converges to zero. Moreover, if
$w^{*}$ is a limit point of sequence $\left\\{w_{k}\right\\}$, then $w^{*}$
satisfies Riemannian KKT conditions (4).

The proof of the above theorem will be given in the next section. Note that
although the exponential map is used in the proof, the numerical experiments
indicate that global convergence may hold for a general retraction
$\mathrm{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
$\mathcal{M}=\mathbb{R}^{n}$ in (RCOP), and we call the algorithm in [11] the
Euclidean Interior Point Method (EIPM). In what follows, we will omit similar
content because of space limitations and focus on the difficulties encountered
when adapting the proof of EIPM to RIPM. In particular, we will make these
difficulties as tractable as in EIPM by proving a series 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 of some special scalar fields on manifold $\mathcal{M}$. The claims
of this subsection are trivial if $\mathcal{M}=\mathbb{R}^{n}$, but they need
to be treated carefully for general $\mathcal{M}$.

If we assign a linear operator $\mathcal{A}_{x}\colon T_{x}\mathcal{M}\to
T_{x}\mathcal{M}$ to each $x\in\mathcal{M},$ then the map
$x\mapsto\left\|\mathcal{A}_{x}\right\|:=\sup\left\\{\left\|\mathcal{A}_{x}v\right\|_{x}\mid
v\in T_{x}\mathcal{M},\|v\|_{x}=1,\text{ or, }\|v\|_{x}\leq 1\right\\}$ is a
scalar field on $\mathcal{M}$; but notice that the operator norm $\|\cdot\|$
depends on $x$. Let $\operatorname{Sym}(d)$ denote the set of symmetric
matrices of order $d$, and $\|\cdot\|_{\mathrm{F}}$, $\|\cdot\|_{2}$ denote
the Frobenius norm and the spectral norm, respectively, applied to a given
matrix.

######  Lemma 7.1.

Let $\mathcal{M}$ be an $n$-dimensional Riemannian manifold. Let
$x\in\mathcal{M}$ and $\mathcal{A}_{x}$ be a linear operator on
$T_{x}\mathcal{M}$. Choose an orthonormal basis of $T_{x}\mathcal{M}$ with
respect to $\langle\cdot,\cdot\rangle_{x}$, and let
$\hat{\mathcal{A}}_{x}\in\mathbb{R}^{n\times n}$ denote the matrix
representation of $\mathcal{A}_{x}$ under the basis. Then,
$\|\hat{\mathcal{A}}_{x}\|_{2}$, $\|\hat{\mathcal{A}}_{x}\|_{\mathrm{F}}$ are
_invariant_ under a change of orthonormal basis; moreover,
$\left\|\mathcal{A}_{x}\right\|=\|\hat{\mathcal{A}}_{x}\|_{2}\leq\|\hat{\mathcal{A}}_{x}\|_{\mathrm{F}}.$

###### Proof.

Suppose that there are two orthonormal bases
$\\{E_{i}\\}_{i=1}^{n},\\{E^{\prime}_{i}\\}_{i=1}^{n}$ on $T_{x}\mathcal{M}$.
With respect to them, let $P\in\mathbb{R}^{n\times n}$ denote the change-of-
basis matrix, i.e., $[P]_{kj}:=\langle E^{\prime}_{j},E_{k}\rangle_{x}$, for
$1\leq k,j\leq n$; then, $P$ is orthogonal. Let
$\hat{\mathcal{A}}_{x},\hat{\mathcal{A}}^{\prime}_{x}\in\mathbb{R}^{n\times
n}$ denote the matrix representations of $\mathcal{A}_{x}$ under the two
bases, respectively. We have
$\hat{\mathcal{A}}_{x}^{\prime}=P^{-1}\hat{\mathcal{A}}_{x}P.$ Then,
$\|\hat{\mathcal{A}}_{x}^{\prime}\|=\|P^{-1}\hat{\mathcal{A}}_{x}P\|=\|\hat{\mathcal{A}}_{x}\|$
holds for the Frobenius norm or the spectral norm. Therefore, the values
$\|\hat{\mathcal{A}}_{x}\|_{2}$ and $\|\hat{\mathcal{A}}_{x}\|_{\mathrm{F}}$
are invariant under a change of orthonormal basis. Now, consider an
orthonormal basis $\\{E_{i}\\}_{i=1}^{n}$ on $T_{x}\mathcal{M}$. For any $y\in
T_{x}\mathcal{M}$, its vector representation $\hat{y}\in\mathbb{R}^{n}$ is
defined by $y=\sum_{i=1}^{n}\hat{y}_{i}E_{i}$. Accordingly, we have
$\widehat{\mathcal{A}_{x}y}=\hat{\mathcal{A}}_{x}\hat{y}$, i.e.,
$\mathcal{A}_{x}y=\sum_{i=1}^{n}(\hat{\mathcal{A}}_{x}\hat{y})_{i}E_{i}$, and
from the orthonormal property of the basis, we have

$\left\|\mathcal{A}_{x}y\right\|_{x}^{2}=\|\sum_{i=1}^{n}(\hat{\mathcal{A}}_{x}\hat{y})_{i}E_{i}\|_{x}^{2}=\sum_{i=1}^{n}(\hat{\mathcal{A}}_{x}\hat{y})_{i}^{2}=\|\hat{\mathcal{A}}_{x}\hat{y}\|_{2}^{2}.$

Thus,
$\left\|\mathcal{A}_{x}y\right\|_{x}=\|\hat{\mathcal{A}}_{x}\hat{y}\|_{2}$ for
any $y\in T_{x}\mathcal{M}$. Finally, we have

$\left\|\mathcal{A}_{x}\right\|=\sup_{\begin{subarray}{c}y\in
T_{x}\mathcal{M},\|y\|_{x}=1\end{subarray}}\left\|\mathcal{A}_{x}y\right\|_{x}=\sup_{\begin{subarray}{c}\hat{y}\in\mathbb{R}^{n},\|\hat{y}\|_{2}=1\end{subarray}}\|\hat{\mathcal{A}}_{x}\hat{y}\|_{2}=\|\hat{\mathcal{A}}_{x}\|_{2}.$

It is clear that $\|X\|_{2}\leq\|X\|_{\mathrm{F}}$ for any matrix $X$. ∎

Yet, we have not clarified the continuity about
$x\mapsto\left\|A_{x}\right\|$. The following proposition proves the
continuity of an important case of $A_{x}$ that appears in (RCOP).

######  Proposition 7.2.

Consider $f$ in (RCOP). Let
$\hat{\operatorname{Hess}}f(x)\in\operatorname{Sym}(d)$ denote the matrix
representation of $\operatorname{Hess}f(x)$ under an arbitrary orthonormal
basis of $T_{x}\mathcal{M}$. Then, $x\mapsto\|\hat{\operatorname{Hess}}f(x)\|$
is a continuous scalar field on $\mathcal{M}$ for $\|\cdot\|_{\mathrm{F}}$ or
$\|\cdot\|_{2}$. Moreover, $x\mapsto\|\operatorname{Hess}f(x)\|$ is a
continuous scalar field on $\mathcal{M}$.

###### Proof.

Lemma 7.1 shows that $x\mapsto\|\hat{\operatorname{Hess}}f(x)\|$ is well
defined, so it suffices to prove its continuity. From Corollary 13.8 in [21],
for each $\bar{x}\in\mathcal{M}$ there is a smooth, orthonormal local frame
$\\{E_{i}\\}_{i=1}^{d}$ on a neighborhood $\mathcal{U}$ of $\bar{x}$; namely,
$\\{E_{1}(x),\ldots,E_{d}(x)\\}$ forms an orthonormal basis on
$T_{x}\mathcal{M}$ for all $x\in\mathcal{U}.$ Choose such a local frame
$\\{E_{i}\\}_{i=1}^{d}$ around $\bar{x}$; then, the matrix representation of
$\operatorname{Hess}f(x)$ is given by
$[\hat{\operatorname{Hess}}f(x)]_{kj}:=\left\langle\operatorname{Hess}f(x)[E_{j}(x)],E_{k}(x)\right\rangle_{x}=\langle(\nabla_{E_{j}}\operatorname{grad}f)(x),E_{k}(x)\rangle_{x}\text{
for }1\leq k,j\leq d.$ From the smoothness of the Riemannian metric (1), it
follows that $x\mapsto\hat{\operatorname{Hess}}f(x)$ is a continuous function
from $\mathcal{U}\subset\mathcal{M}$ to $\operatorname{Sym}(d)$. Since matrix
norms are continuous, $\|\hat{\operatorname{Hess}}f(x)\|$ is continuous on
$\mathcal{U}\ni\bar{x}.$ This argument holds for any $\bar{x}\in\mathcal{M}$.
From Lemma 7.1,
$\|\hat{\operatorname{Hess}}f(x)\|_{2}=\|\operatorname{Hess}f(x)\|$ for any
$x\in\mathcal{M}$, which completes the proof. ∎

The above result can be applied verbatim to the Hessian of constraint
functions $\left\\{h_{i}\right\\}_{i=1}^{l}$,
$\left\\{g_{i}\right\\}_{i=1}^{m}$ in (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 $\mathcal{H}_{x}$,
$\mathcal{G}_{x}$ defined in (7). Then,
$x\mapsto\left\|\mathcal{H}_{x}\right\|$ and
$x\mapsto\left\|\mathcal{G}_{x}\right\|$ are continuous scalar fields on
$\mathcal{M}$.

######  Proposition 7.4.

Given $w=(x,y,z,s)\in\mathcal{N}$, consider the operator $\nabla F(w)$ in (8).
Let $\\{E_{i}\\}_{i=1}^{d}$ be an orthonormal basis of $T_{x}\mathcal{M}$ and
$\\{e_{i}\\}_{i=1}^{l}$, $\\{\dot{e}_{i}\\}_{i=1}^{m}$ be the standard bases
of $\mathbb{R}^{l},\mathbb{R}^{m}$, respectively. If we choose an orthonormal
basis of $T_{w}\mathcal{N}$ as follows:

$\left\\{\left(E_{i},0,0,0\right)\right\\}_{i=1}^{d}\cup\left\\{\left(0_{x},e_{i},0,0\right)\right\\}_{i=1}^{l}\cup\left\\{\left(0_{x},0,\dot{e}_{i},0\right)\right\\}_{i=1}^{m}\cup\left\\{\left(0_{x},0,0,\dot{e}_{i}\right)\right\\}_{i=1}^{m},$
(35)

then, the matrix representation of $\nabla F(w)$ is given by

$\hat{\nabla}F(w)=\left(\begin{array}[]{cccc}Q&B&C&0\\\ B^{T}&0&0&0\\\
C^{T}&0&0&I\\\ 0&0&S&Z\end{array}\right),$

i.e., a matrix of order $(d+l+2m)$ and, where

$Q:=Q(w)\in\operatorname{Sym}(d)$ is given by
$[Q]_{kj}:=\left\langle\operatorname{Hess}_{x}\mathcal{L}(w)[E_{j}],E_{k}\right\rangle_{x}$
for $1\leq k,j\leq d$;

$B:=B(x)=[\hat{\operatorname{grad}}h_{1}(x),\cdots,\hat{\operatorname{grad}}h_{l}(x)]\in\mathbb{R}^{d\times
l};$

$C:=C(x)=[\hat{\operatorname{grad}}g_{1}(x),\cdots,\hat{\operatorname{grad}}g_{m}(x)]\in\mathbb{R}^{d\times
m}$ and the “hat” symbol above means the corresponding vector representation
under the basis $\\{E_{i}\\}_{i=1}^{d}$.

In this case, there is a continuous scalar field
$T\colon\mathcal{N}\to\mathbb{R}$ such that for any $w$,
$\|Q(w)\|_{\mathrm{F}}\leq T(w)$. Moreover, $x\mapsto\|B(x)\|_{\mathrm{F}}$
and $x\mapsto\|C(x)\|_{\mathrm{F}}$ are continuous scalar fields on
$\mathcal{M}$.

###### Proof.

The matrix $\hat{\nabla}F(w)$ under the basis (35) is obtained through a
trivial process, so we will omit its description. From relation (9), we have
$Q(w):=\hat{\operatorname{Hess}}_{x}\mathcal{L}(w)=\hat{\operatorname{Hess}}f(x)+\sum_{i=1}^{l}y_{i}\hat{\operatorname{Hess}}h_{i}(x)+\sum_{i=1}^{m}z_{i}\hat{\operatorname{Hess}}g_{i}(x),$
under the same basis. Thus,

$\left\|Q(w)\right\|_{\mathrm{F}}\leq\|\hat{\operatorname{Hess}}f(x)\|_{\mathrm{F}}+\sum_{i=1}^{l}|y_{i}|\|\hat{\operatorname{Hess}}h_{i}(x)\|_{\mathrm{F}}+\sum_{i=1}^{m}|z_{i}|\|\hat{\operatorname{Hess}}g_{i}(x)\|_{\mathrm{F}}=:T(w).$

From Proposition 7.2, $\|\hat{\operatorname{Hess}}f(x)\|_{\mathrm{F}}$,
$\|\hat{\operatorname{Hess}}h_{i}(x)\|_{\mathrm{F}}$, and
$\|\hat{\operatorname{Hess}}g_{i}(x)\|_{\mathrm{F}}$ are all continuous with
respect to $x$, thus, $T$ is continuous. As for
$x\mapsto\|B(x)\|_{\mathrm{F}}$, since the basis $\\{E_{i}\\}_{i=1}^{d}$ is
orthonormal,

$\|B(x)\|_{\mathrm{F}}^{2}=\sum_{i=1}^{l}\|\hat{\operatorname{grad}}h_{i}(x)\|_{2}^{2}=\sum_{i=1}^{l}\left\|\operatorname{grad}h_{i}(x)\right\|_{x}^{2},$

which implies continuity by (1). The claim for $x\mapsto\|C(x)\|_{\mathrm{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 difficulties encountered when adapting the proof of EIPM to RIPM.

######  Proposition 7.5 (Boundedness of the sequences).

Let $\left\\{w_{k}\right\\}$ be a sequence generated by Algorithm 5 and
suppose that (C1)-(C3) hold. If $\epsilon>0$ and $w_{k}\in\Omega(\epsilon)$
for all $k$, then

  1. (a)

$\\{z_{k}^{T}s_{k}\\}$,
$\left\\{(z_{k})_{i}(s_{k})_{i}\right\\},i=1,\ldots,m$, are all bounded above
and below away from zero;

  2. (b)

$\left\\{z_{k}\right\\}$ and $\left\\{s_{k}\right\\}$ are bounded above and
component-wise bounded away from zero;

  3. (c)

$\left\\{w_{k}\right\\}$ is bounded;

  4. (d)

$\\{\|\nabla F(w_{k})^{-1}\|\\}$ is bounded;

  5. (e)

$\left\\{\Delta w_{k}\right\\}$ 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 suffices to prove that $\left\\{y_{k}\right\\}$ is bounded. The
iteration count $k$ is omitted in what follows. By using the notation
$\mathcal{H}_{x}$ and $\mathcal{G}_{x}$ as defined in (7), we have
$\mathcal{H}_{x}y=\operatorname{grad}_{x}\mathcal{L}(w)-\operatorname{grad}f(x)-\mathcal{G}_{x}z=:b.$
By (C1), $\mathcal{H}_{x}$ is an injection; then, there exists a unique
solution $y$ to $\mathcal{H}_{x}y=b$. Indeed, we have

$y=[\left(\mathcal{H}_{x}^{*}\mathcal{H}_{x}\right)^{-1}\mathcal{H}_{x}^{*}]\left(\operatorname{grad}_{x}\mathcal{L}(w)-\operatorname{grad}f(x)-\mathcal{G}_{x}z\right).$
(36)

Define $\mathcal{C}_{x}:T_{x}\mathcal{M}\to\mathbb{R}^{l}$ as
$\mathcal{C}_{x}:=\left(\mathcal{H}_{x}^{*}\mathcal{H}_{x}\right)^{-1}\mathcal{H}_{x}^{*}.$
Under an arbitrary orthonormal basis of $T_{x}\mathcal{M}$ and standard basis
of $\mathbb{R}^{l}$, if $\hat{\mathcal{H}}_{x}$ is the matrix corresponding to
$\mathcal{H}_{x}$, then
$\hat{\mathcal{C}}_{x}=(\hat{\mathcal{H}}_{x}^{T}\hat{\mathcal{H}}_{x})^{-1}\hat{\mathcal{H}}_{x}^{T}$.
It is easy to show that
$\left\|\mathcal{C}_{x}\right\|=\|\hat{\mathcal{C}}_{x}\|_{2}$ for any $x$. We
can see that for each $\bar{x}\in\mathcal{M}$ there is a neighborhood
$\mathcal{U}$ of $\bar{x}$ such that $x\mapsto\hat{\mathcal{H}}_{x}$ is
continuous over $\mathcal{U}$. Then, by function composition,
$x\mapsto\hat{\mathcal{C}}_{x}$ is also continuous over $\mathcal{U}$. This
implies that
$x\mapsto\left\|\mathcal{C}_{x}\right\|=\|\hat{\mathcal{C}}_{x}\|_{2}$ is
continuous at each $\bar{x}$, and hence, on $\mathcal{M}$. Finally, by
Proposition 7.3,
$\left\|\mathcal{C}_{x}\right\|,\left\|\operatorname{grad}f(x)\right\|,$ and
$\left\|\mathcal{G}_{x}\right\|$ are all continuous on $\mathcal{M}$. Because
$\\{x_{k}\\}$ is bounded, by (36) we have
$\left\|y_{k}\right\|\leq\left\|\mathcal{C}_{x_{k}}\right\|\left(\left\|\operatorname{grad}_{x}\mathcal{L}(w_{k})\right\|+\left\|\operatorname{grad}f(x_{k})\right\|+\left\|\mathcal{G}_{x_{k}}\right\|\left\|z_{k}\right\|\right)\leq
c_{1}\left(\sqrt{\varphi_{0}}+c_{2}+c_{3}\left\|z_{k}\right\|\right),$ for
some positive constants $c_{1},c_{2},c_{3}.$ Then, $\left\\{y_{k}\right\\}$ is
bounded because $\left\\{z_{k}\right\\}$ is bounded.

(d) For each $w_{k}$, choose an arbitrary orthonormal basis of
$T_{w_{k}}\mathcal{N}.$ If the matrix representation $\hat{\nabla}F(w_{k})$
corresponds to $\nabla F(w_{k})$, then $[\hat{\nabla}F(w_{k})]^{-1}$
corresponds to $\nabla F(w_{k})^{-1}$. By Lemma 7.1, we have $\left\|\nabla
F(w_{k})^{-1}\right\|\leq\|[\hat{\nabla}F(w_{k})]^{-1}\|_{\mathrm{F}};$ thus,
it is sufficient to show that
$\\{\|[\hat{\nabla}F(w_{k})]^{-1}\|_{\mathrm{F}}\\}$ is bounded. For
convenience, we will choose the basis of $T_{w_{k}}\mathcal{N}$ given in (35).
Then, we have

$\hat{\nabla}F(w_{k})=\left(\begin{array}[]{cccc}Q_{k}&B_{k}&C_{k}&0\\\
B^{T}_{k}&0&0&0\\\ C^{T}_{k}&0&0&I\\\ 0&0&S_{k}&Z_{k}\end{array}\right).$

By Proposition 7.4, there is a continuous scalar field
$T:\mathcal{N}\rightarrow\mathbb{R}$ such that $\|Q(w)\|_{\mathrm{F}}\leq
T(w)$ for all $w\in\mathcal{N}$; and $\|B(x)\|_{\mathrm{F}}$,
$\|C(x)\|_{\mathrm{F}}$ are continuous on $\mathcal{M}$. It follows from the
boundedness of $\left\\{x_{k}\right\\}$ and $\left\\{w_{k}\right\\}$ that for
all $k$,
$\left\|Q_{k}\right\|_{\mathrm{F}}\equiv\left\|Q(w_{k})\right\|_{\mathrm{F}}\leq
T(w_{k})\leq c_{4},$
$\left\|B_{k}\right\|_{\mathrm{F}}\equiv\left\|B(x_{k})\right\|_{\mathrm{F}}\leq
c_{5},$ and
$\left\|C_{k}\right\|_{\mathrm{F}}\equiv\left\|C(x_{k})\right\|_{\mathrm{F}}\leq
c_{6},$ for some positive constants $c_{4},c_{5},$ and $c_{6}.$

On the other hand, whichever basis is used in the form of (35), the structure
of $\hat{\nabla}F(w_{k})$ and the properties of its block submatrices remain
unchanged, e.g., symmetry of $Q_{k}$; full rank of $B_{k}$; identity matrix
$I$ in the third row; all zero matrices; diagonal matrices $S_{k},Z_{k}$; etc.
This ensures that we can obtain the desired result by performing an
appropriate decomposition of $\hat{\nabla}F(w_{k})$. Up to this point, we have
created all the conditions needed in the proof of the Euclidean version,
namely, EIPM. We can make the claim that
$\\{\|[\hat{\nabla}F(w_{k})]^{-1}\|_{\mathrm{F}}\\}$ is bounded by applying
the proofs in [11, Lemma 6.2] and/or [3, Theorem 2 (c)] directly. ∎

######  Lemma 7.6 ($\left\\{\bar{\alpha}_{k}\right\\}$ bounded away from
zero).

Let $\left\\{w_{k}\right\\}$ be generated by Algorithm 5 with
$\mathrm{R}=\operatorname{Exp}$ and let (C1)-(C3) hold. If $\epsilon>0$ and
$w_{k}\in\Omega(\epsilon)$ for all $k$, $\left\\{\sigma_{k}\right\\}$ is
bounded away from zero; then, $\left\\{\bar{\alpha}_{k}\right\\}$ is bounded
away from zero.

###### Proof.

Since $\bar{\alpha}_{k}=\min\\{\alpha_{k}^{I},\alpha_{k}^{II}\\}$, it is
sufficient to show that $\\{\alpha_{k}^{I}\\}$ and $\\{\alpha_{k}^{II}\\}$ are
bounded away from zero. For $\alpha_{k}^{I}$, 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 $\alpha_{k}^{II}$, we need to adapt the proofs in
[11, 9], since Lipschitz continuity on manifolds is more complicated, see
Subsection 4.2.

Let us suppress the subscript $k$. Recall that
$w(\alpha)=\bar{\operatorname{Exp}}_{w}(\alpha\Delta w).$ Fix $\alpha\Delta w$
and let $\mathrm{P}_{\gamma}$ be the parallel transport along the geodesic
$\gamma(t)=\bar{\operatorname{Exp}}_{w}\left(t\alpha\Delta w\right)$. By Lemma
4.8 where $c_{1}=0$, we obtain

$\displaystyle\mathrm{P}_{\gamma}^{0\to 1}[F(w(\alpha))]$
$\displaystyle=F(w)+\alpha\nabla F(w)[\Delta w]-\alpha\nabla F(w)[\Delta
w]+\int_{0}^{1}\mathrm{P}_{\gamma}^{0\to t}\nabla
F(\gamma(t))\mathrm{P}_{\gamma}^{t\to 0}[\alpha\Delta w]\mathrm{d}t$
$\displaystyle=F(w)+\alpha\left(\sigma\rho\hat{e}-F(w)\right)+\alpha\int_{0}^{1}\left(\mathrm{P}_{\gamma}^{0\to
t}\nabla F(\gamma(t))\mathrm{P}_{\gamma}^{t\to 0}-\nabla F(w)\right)[\Delta
w]\mathrm{d}t$
$\displaystyle=(1-\alpha)F(w)+\alpha\sigma\rho\hat{e}+\alpha\int_{0}^{1}\left(\mathrm{P}_{\gamma}^{0\to
t}\nabla F(\gamma(t))\mathrm{P}_{\gamma}^{t\to 0}-\nabla F(w)\right)[\Delta
w]\mathrm{d}t.$

Taking the norm on both sides above gives

$\displaystyle\left\|F(w(\alpha))\right\|$
$\displaystyle=\left\|\mathrm{P}_{\gamma}^{0\to 1}[F(w(\alpha))]\right\|\text{
(since parallel transport (\ref{eq:ParallelTransport}) is isometric)}$
$\displaystyle\leq(1-\alpha)\left\|F(w)\right\|+\alpha\sigma\rho\left\|\hat{e}\right\|+\alpha\int_{0}^{1}\left\|\mathrm{P}_{\gamma}^{0\to
t}\nabla F(\gamma(t))\mathrm{P}_{\gamma}^{t\to 0}-\nabla
F(w)\right\|\left\|\Delta w\right\|\mathrm{d}t$
$\displaystyle\leq(1-\alpha)\left\|F(w)\right\|+\alpha\sigma\rho\sqrt{m}+\alpha\int_{0}^{1}\kappa_{2}\|t\alpha\Delta
w\|\left\|\Delta w\right\|\mathrm{d}t$
$\displaystyle=(1-\alpha)\left\|F(w)\right\|+\alpha\sigma\rho\sqrt{m}+\alpha^{2}\|\Delta
w\|^{2}\kappa/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 $\left\\{\left\|F(w_{k})\right\|\right\\}$ is
monotonically nonincreasing, hence convergent. Assume that
$\left\\{\left\|F(w_{k})\right\|\right\\}$ does not converge to zero. Then,
there exists $\epsilon>0$ such that
$\left\\{w_{k}\right\\}\subset\Omega(\epsilon)$ for infinitely many $k$. We
will show that the following two cases both lead to contradictions and thus,
the hypothesis $\left\|F(w_{k})\right\|\nrightarrow 0$ is not valid.

Case 1. For infinitely many $k$, if step (3b) in Algorithm 5 is executed with
$\alpha_{k}\equiv\bar{\alpha}_{k}$, it follows from Proposition 6.2 that
$\varphi(w_{k+1})/\varphi(w_{k})\leq\lambda_{k}:=\left[1-2\bar{\alpha}_{k}\beta\left(1-\sigma_{k}\right)\right].$
Since $\left\\{\bar{\alpha}_{k}\right\\}$ is bounded away from zero by Lemma
7.6 and $\left\\{\sigma_{k}\right\\}$ is bounded away from one, then
$\left\\{\lambda_{k}\right\\}$ is bounded away from one and hence,
$\varphi(w_{k})\to 0$; this is a contradiction.

Case 2. On the other hand, for infinitely many $k$, if
$\alpha_{k}<\bar{\alpha}_{k}$, we have that
$\alpha_{k}\leq\theta\bar{\alpha}_{k}$. Then, condition (33) fails to hold for
an $\tilde{\alpha}_{k}$ with
$\alpha_{k}<\tilde{\alpha}_{k}\leq\alpha_{k}/\theta=\theta^{t-1}\bar{\alpha}_{k}$.
Notice that $\alpha_{k}/\theta$ is the value corresponding to the last
failure. Recall that the derivative of the real-valued function
$\alpha\mapsto\varphi(\alpha):=\varphi\left(\bar{\operatorname{Exp}}_{w_{k}}(\alpha\Delta
w_{k})\right)$ is

$\varphi^{\prime}(\alpha)=\mathcal{D}\varphi(\bar{\operatorname{Exp}}_{w_{k}}\left(\alpha\Delta
w_{k}\right))[\mathcal{D}\bar{\operatorname{Exp}}_{w_{k}}\left(\alpha\Delta
w_{k}\right)\left[\Delta w_{k}\right]].$ (37)

Applying the mean value theorem to $\varphi(\alpha)$ on the interval
$[0,\tilde{\alpha}_{k}]$ yields a number $\xi\in(0,1)$ such that
$\tilde{\alpha}_{k}\varphi^{\prime}(\xi\tilde{\alpha}_{k})=\varphi(\tilde{\alpha}_{k})-\varphi(0).$
For short, let $u:=\xi\tilde{\alpha}_{k}\Delta w_{k}.$ Hence,

$\displaystyle\tilde{\alpha}_{k}\beta\left\langle\operatorname{grad}\varphi_{k},\Delta
w_{k}\right\rangle$ (38)
$\displaystyle<\varphi(\tilde{\alpha}_{k})-\varphi(0)\text{ (as condition
(\ref{eq:ArmijoCondition}) fails for $\tilde{\alpha}_{k}$)}$
$\displaystyle=\tilde{\alpha}_{k}\varphi^{\prime}(\xi\tilde{\alpha}_{k})$
$\displaystyle=\tilde{\alpha}_{k}\mathcal{D}\varphi(\bar{\operatorname{Exp}}_{w_{k}}(u))\left[\mathcal{D}\bar{\operatorname{Exp}}_{w_{k}}(u)\left[\Delta
w_{k}\right]\right]\text{ (by equation (\ref{eq:241}))}$
$\displaystyle=\tilde{\alpha}_{k}\langle\operatorname{grad}\varphi(\bar{\operatorname{Exp}}_{w_{k}}(u)),\mathcal{D}\bar{\operatorname{Exp}}_{w_{k}}(u)\left[\Delta
w_{k}\right]\rangle.$

On the other hand, note that

$\displaystyle\left\langle\operatorname{grad}\varphi_{k},\Delta
w_{k}\right\rangle$ (39)
$\displaystyle=\langle\operatorname{grad}\varphi_{k},u\rangle/\xi\tilde{\alpha}_{k}$
$\displaystyle=\langle\mathcal{D}\bar{\operatorname{Exp}}_{w_{k}}(u)\left[\operatorname{grad}\varphi_{k}\right],\mathcal{D}\bar{\operatorname{Exp}}_{w_{k}}(u)\left[u\right]\rangle/\xi\tilde{\alpha}_{k}\text{
(by Lemma \ref{lem:gauss})}$
$\displaystyle=\langle\mathcal{D}\bar{\operatorname{Exp}}_{w_{k}}(u)\left[\operatorname{grad}\varphi_{k}\right],\mathcal{D}\bar{\operatorname{Exp}}_{w_{k}}(u)\left[\Delta
w_{k}\right]\rangle.$

Subtracting
$\tilde{\alpha}_{k}\left\langle\operatorname{grad}\varphi_{k},\Delta
w_{k}\right\rangle$ from both sides of (38) and using equalities (39) gives

$\displaystyle\tilde{\alpha}_{k}(\beta-1)\left\langle\operatorname{grad}\varphi_{k},\Delta
w_{k}\right\rangle$
$\displaystyle<\tilde{\alpha}_{k}\left[\langle\operatorname{grad}\varphi(\bar{\operatorname{Exp}}_{w_{k}}(u)),\mathcal{D}\bar{\operatorname{Exp}}_{w_{k}}(u)\left[\Delta
w_{k}\right]\rangle-\langle\operatorname{grad}\varphi_{k},\Delta
w_{k}\rangle\right]$
$\displaystyle=\tilde{\alpha}_{k}\langle\operatorname{grad}\varphi(\bar{\operatorname{Exp}}_{w_{k}}(u))-\mathcal{D}\bar{\operatorname{Exp}}_{w_{k}}(u)\left[\operatorname{grad}\varphi_{k}\right],\mathcal{D}\bar{\operatorname{Exp}}_{w_{k}}(u)\left[\Delta
w_{k}\right]\rangle$
$\displaystyle\leq\tilde{\alpha}_{k}\left\|\operatorname{grad}\varphi(\bar{\operatorname{Exp}}_{w_{k}}(u))-\mathcal{D}\bar{\operatorname{Exp}}_{w_{k}}(u)\left[\operatorname{grad}\varphi_{k}\right]\right\|\left\|\mathcal{D}\bar{\operatorname{Exp}}_{w_{k}}(u)\left[\Delta
w_{k}\right]\right\|$
$\displaystyle=\tilde{\alpha}_{k}\left\|\operatorname{grad}\varphi(y)-\mathcal{D}\bar{\operatorname{Exp}}_{w_{k}}(u)\left[\operatorname{grad}\varphi(w_{k})\right]\right\|\left\|\mathcal{D}\bar{\operatorname{Exp}}_{w_{k}}(u)\left[\Delta
w_{k}\right]\right\|$ (by letting $y:=\bar{\operatorname{Exp}}_{w_{k}}(u)$)
$\displaystyle\leq\tilde{\alpha}_{k}\kappa\left\|u\right\|\left\|\Delta
w_{k}\right\|\text{ (by (\ref{eq:DifferentiatedRetraction}), Definition
\ref{defn:LipschitzContinuouslyDifferentiable} and Lemma \ref{lem:gauss})}$
$\displaystyle=\kappa\xi\tilde{\alpha}_{k}^{2}\left\|\Delta
w_{k}\right\|^{2}.$

Finally, we obtain $(\beta-1)\left\langle\operatorname{grad}\varphi_{k},\Delta
w_{k}\right\rangle/(\kappa\xi\left\|\Delta
w_{k}\right\|^{2})<\tilde{\alpha}_{k}.$

Because $\left\langle\operatorname{grad}\varphi_{k},\Delta
w_{k}\right\rangle<0$ and $\alpha_{k}$ satisfies (33), we have

$\displaystyle\varphi_{k}(0)-\varphi_{k}(\alpha_{k})$
$\displaystyle\geq-\alpha_{k}\beta\left\langle\operatorname{grad}\varphi_{k},\Delta
w_{k}\right\rangle$
$\displaystyle\geq-\theta\beta\tilde{\alpha}_{k}\left\langle\operatorname{grad}\varphi_{k},\Delta
w_{k}\right\rangle$
$\displaystyle\geq-\theta\beta\left\langle\operatorname{grad}\varphi_{k},\Delta
w_{k}\right\rangle(\beta-1)\left\langle\operatorname{grad}\varphi_{k},\Delta
w_{k}\right\rangle/(\kappa\xi\left\|\Delta w_{k}\right\|^{2})$
$\displaystyle\geq[\theta\beta(1-\beta)/\kappa\xi]\left(\left\langle\operatorname{grad}\varphi_{k},\Delta
w_{k}\right\rangle/\left\|\Delta w_{k}\right\|\right)^{2}$
$\displaystyle=\omega\left(\left\langle\operatorname{grad}\varphi_{k},\Delta
w_{k}\right\rangle/\left\|\Delta w_{k}\right\|\right),$

where $\omega(\cdot)$ is an $F$-function (see [28, Definition 14.2.1 & 14.2.2
in P479]). Since $\\{\varphi_{k}\\}$ is bounded below and
$\varphi_{k}\geq\varphi_{k+1}$, it follows that
$\lim_{k\to\infty}(\varphi_{k}-\varphi_{k+1})=0$. By the definition of
$F$-functions, we obtain $\left\langle\operatorname{grad}\varphi_{k},\Delta
w_{k}\right\rangle/\left\|\Delta w_{k}\right\|\to 0.$ Since $\\{\left\|\Delta
w_{k}\right\|\\}$ is bounded (Proposition 7.5), we have
$\left\langle\operatorname{grad}\varphi_{k},\Delta w_{k}\right\rangle\to 0.$
Choosing $\rho_{k}$ with
$z_{k}^{T}s_{k}/m\leq\rho_{k}\leq\left\|F(w_{k})\right\|/\sqrt{m}$ in
Algorithm 5 implies that

$\displaystyle\left\langle\operatorname{grad}\varphi_{k},\Delta
w_{k}\right\rangle/(-2)$
$\displaystyle=\varphi_{k}-\sigma_{k}\rho_{k}z_{k}^{T}s_{k}\geq\varphi_{k}-\sigma_{k}\left\|F(w_{k})\right\|z_{k}^{T}s_{k}/\sqrt{m}$
$\displaystyle\geq\varphi_{k}-\sigma_{k}\left\|F(w_{k})\right\|^{2}\geq(1-\sigma_{k})\varphi_{k}.$

This shows that $\varphi(w_{k})\to 0$, because $\left\\{\sigma_{k}\right\\}$
is bounded away from one; this is a contradiction. ∎

##  8 Numerical Experiments

The numerical experiments compared the performance of the globally convergent
RIPM (Algorithm 5) 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:

  * •

fixed-rank manifold, $\mathcal{M}_{r}=\left\\{X\in\mathbb{R}^{m\times
n}:\operatorname{rank}(X)=r\right\\}$;

  * •

Stiefel manifold, $\operatorname{St}(n,k)=\left\\{X\in\mathbb{R}^{n\times
k}:X^{T}X=I_{k}\right\\}$;

  * •

oblique manifold, $\mathrm{Ob}(n,k)=\\{X\in\mathbb{R}^{n\times
k}:\left(X^{T}X\right)_{ii}=1,\forall i=1,\dots,k\\}.$

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 solve

$\min_{X\in\mathcal{M}_{r}}\|A-X\|_{\mathrm{F}}^{2}\quad\text{ s.t. }X\geq 0.$
(NLRM)

_Input._ We tested three cases of integer triples $(m,n,r)$: $(20,16,2)$,
$(30,24,3)$, and $(40,32,4)$. For each $(m,n,r)$, we generated nonnegative
$L\in\mathbb{R}^{m\times r}$, $R\in\mathbb{R}^{r\times n}$ whose entries
follow a uniform distribution in [0,1]. Then, we added the Gaussian noise with
zero mean and different standard deviation $(\sigma=0,0.001,0.01)$ to $A=LR$.
When there is no noise (i.e., $\sigma=0$), the input data matrix $A$ itself is
exactly a solution.

Problem II. Given $C\in\mathbb{R}^{n\times k}$, [19] computed its projection
onto the nonnegative part of the Stiefel manifold. If the distance is measured
in terms of $\|C-X\|_{\mathrm{F}}^{2}$, we can express it equivalently as

$\min_{X\in\operatorname{St}(n,k)}-2\operatorname{trace}(X^{T}C)\quad\text{
s.t. }X\geq 0.$ (Model_St)

In [19], it is shown that (Model_St) can be equivalently reformulated into

$\min_{X\in\operatorname{Ob}(n,k)}-2\operatorname{trace}(X^{T}C)\quad\text{
s.t. }X\geq 0,\;\|XV\|_{\mathrm{F}}=1,$ (Model_Ob)

where the positive integer $p$ and $V\in\mathbb{R}^{k\times p}$ can be
arbitrary as long as $\|V\|_{\mathrm{F}}=1$ and $VV^{T}$ is entry-wise
positive. We examined both models. _Input._ We considered four cases of
integer pairs $(n,k)$: $(40,8),(50,10),(60,16),$ and $(70,14).$ For a general
$C$, it is always difficult to seek nonnegative projections globally.
Fortunately, Proposition 1 in [19] showed a way to construct $C$ such that
Model_St has a 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 $z_{0}$
and $s_{0}$ from a uniform distribution in [0,1] and set $y_{0}=0$ if equality
constraints exist. We used $\rho_{k}=z_{k}^{T}s_{k}/m$,
$\sigma_{k}=\min\\{0.5,\|F(w_{k})\|^{1/2}\\}$ 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
sufficient decreasing conditions. We set $\gamma_{-1}=0.9$,
$\gamma_{k+1}=(\gamma_{k}+0.5)/2;$ and $\beta=10^{-4}$, $\theta=0.5$. We
compared RIPM with the following Riemannian methods:

  * •

RALM: Riemannian augmented Lagrangian method [23].

  * •

REPM${}_{\text{lqh}}$: Riemannian exact penalty method with smoothing
functions of linear-quadratic and pseudo-Huber [23].

  * •

REPM${}_{\text{lse}}$: Riemannian exact penalty method with smoothing
functions of log-sum-exp [23].

  * •

RSQP: Riemannian sequential quadratic programming [27].

Let $[t]_{+}:=\max(0,t)$ and $[t]_{-}:=\min(0,t)$ for any $t\in\mathbb{R}.$
The experimental settings followed those of [27], where they used the KKT
residual defined as

$\sqrt{\left\|\operatorname{grad}_{x}\mathcal{L}(w)\right\|^{2}+\sum_{i=1}^{m}\\{[z_{i}]_{-}^{2}+[g_{i}(x)]_{+}^{2}+\left|z_{i}g_{i}(x)\right|^{2}\\}+\sum_{i=1}^{l}\left|h_{i}(x)\right|^{2}}$

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 the same initial point. The experiment is considered successfully
terminated if it finds a solution with a KKT residual lower than
$\epsilon_{\text{kkt}}$ before triggering any of the stopping conditions. For
the _first-order_ algorithms (including RALM and the REPMs), the stopping
conditions are: elapsed time exceeding $t_{\max}$ 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 $t_{\max}$ 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 sufficiently large value. The selection of
$t_{\max}$ is related to the actual time it took to run all the codes on the
computer. Setting $t_{\max}$ too large resulted in excessive time spent on
poorly performing algorithms. On the other hand, $\epsilon_{\mathrm{kkt}}$ was
chosen to better demonstrate that second-order algorithms could achieve more
accurate solutions. Therefore, we chose the appropriate values for $t_{\max}$
and $\epsilon_{\text{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 the average iteration number (Iter.) among the successful
trials. In order to capture the combination of stability and speed, the
algorithms with 0.9 or higher are first 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 same
as the “theoretical convergence in any global convergence theorem”. Numerical
experiments reflect the actual performance of the algorithm in the
application. For example, the first-order algorithms (RALM, REPMs) also have
theoretical global convergence, but under our high criterion, it is difficult
for them to generate a high-precision solution within a certain period of
time. The second-order algorithm (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 $t_{\max}=180$, $\epsilon_{\text{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 success
rates dropped sharply as the noise level (standard deviation $\sigma$)
intensified, eventually leading to non-convergence. In contrast, RSQP and RIPM
were more stable and robust, while RIPM was much faster than RSQP. The cost of
RSQP increased drastically with the problem size because a _quadratic
programming subproblem_ (defined over tangent space of current point) had to
be solved in each iteration. Unlike RIPM, which uses Krylov subspace methods
introduced in Subsection 3.4.2 to avoid 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 same
amount of time as RALM and the REPMs did.

Table 1: Performance of various Riemannian methods on problem I (NLRM). $(m,n,r)$ | (20,16,2) | (30,24,3) | (40,32,4)
---|---|---|---
no noise | Success | Time [s] | Iter. | Success | Time [s] | Iter. | Success | Time [s] | Iter.
RALM | 0.4 | 1.115 | 31 | 0.65 | 1.813 | 31 | 0.75 | 2.800 | 31
REPM${}_{\text{lqh}}$ | 1 | 5.165$\times 10^{-1}$ | 31 | 1 | 1.009 | 31 | 1 | 1.747 | 31
REPM${}_{\text{lse}}$ | 1 | 2.242 | 31 | 1 | 4.041 | 31 | 0.95 | 6.952 | 31
RSQP | 0.9 | 6.429 | 7 | 0.9 | 3.944$\times 10$ | 8 | 0.9 | 1.254$\times 10^{2}$ | 8
RIPM | 1 | $\mathbf{4.920\times 10^{-1}}$ | 19 | 1 | 2.247 | 27 | 1 | 5.277 | 32
$(m,n,r)$ | (20,16,2) | (30,24,3) | (40,32,4)
$\sigma=0.001$ | Success | Time [s] | Iter. | Success | Time [s] | Iter. | Success | Time [s] | Iter.
RALM | 0.2 | 1.001 | 31 | 0.15 | 2.050 | 31 | 0.05 | 2.758 | 31
REPM${}_{\text{lqh}}$ | 0.1 | 4.983$\times 10^{-1}$ | 32 | 0.25 | 1.035 | 31 | 0.15 | 1.787 | 31
REPM${}_{\text{lse}}$ | 0.15 | 2.444 | 31 | 0.1 | 4.867 | 31 | 0.05 | 8.371 | 31
RSQP | 0.95 | 6.619 | 7 | 0.95 | 3.848$\times 10$ | 8 | 0.9 | 1.299$\times 10^{2}$ | 8
RIPM | 1 | $\mathbf{5.376\times 10^{-1}}$ | 20 | 1 | 2.342 | 27 | 1 | 4.631 | 29
$(m,n,r)$ | (20,16,2) | (30,24,3) | (40,32,4)
$\sigma=0.01$ | Success | Time [s] | Iter. | Success | Time [s] | Iter. | Success | Time [s] | Iter.
RALM | 0 | - | - | 0 | - | - | 0 | - | -
REPM${}_{\text{lqh}}$ | 0 | - | - | 0 | - | - | 0 | - | -
REPM${}_{\text{lse}}$ | 0 | - | - | 0 | - | - | 0 | - | -
RSQP | 1 | 7.295 | 8 | 0.95 | 4.114$\times 10$ | 8 | 0.95 | 1.430$\times 10^{2}$ | 9
RIPM | 1 | 5.980$\times 10^{-1}$ | 21 | 0.95 | 1.883 | 25 | 0.95 | 4.602 | 29

Problem II. Here, we set $t_{\max}=600$, $\epsilon_{\text{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 $\|\tilde{X}-X^{*}\|_{\mathrm{F}}$,
where $\tilde{X}$ denotes the final iterate. The numerical results are listed
in Table 2 and 3. The Error columns show that if the KKT residual is
sufficiently small, then $\tilde{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, the success 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 in some cases. Overall, our RIPM was relatively fast and most
stable.

Table 2: Performance of various Riemannian methods on problem II (Model_St). $(n,k)$ | (40,8) | (50,10)
---|---|---
| Success | Time [s] | Iter. | Error | Success | Time [s] | Iter. | Error
RALM | 1 | 2.347 | 45 | 5.41$\times 10^{-7}$ | 1 | 4.344 | 54 | 5.21$\times 10^{-7}$
REPM${}_{\text{lqh}}$ | 0 | - | - | - | 0 | - | - | -
REPM${}_{\text{lse}}$ | 0 | - | - | - | 0 | - | - | -
RSQP | 0.9 | 1.352$\times 10$ | 7 | 2.05$\times 10^{-9}$ | 0.7 | 3.097$\times 10$ | 6 | 2.47$\times 10^{-9}$
RIPM | 1 | 2.225 | 31 | 3.72$\times 10^{-8}$ | 1 | 3.785 | 32 | 3.38$\times 10^{-8}$
$(n,k)$ | (60,12) | (70,14)
| Success | Time [s] | Iter. | Error | Success | Time [s] | Iter. | Error
RALM | 1 | 4.097 | 34 | 4.93$\times 10^{-7}$ | 1 | 6.234 | 37 | 5.34$\times 10^{-7}$
REPM${}_{\text{lqh}}$ | 0 | - | - | - | 0 | - | - | -
REPM${}_{\text{lse}}$ | 0 | - | - | - | 0 | - | - | -
RSQP | 0.65 | 7.802$\times 10$ | 7 | 6.48$\times 10^{-9}$ | 0.85 | 1.661$\times 10^{2}$ | 7 | 2.64$\times 10^{-9}$
RIPM | 1 | 5.555 | 32 | 2.81$\times 10^{-8}$ | 1 | 7.574 | 33 | 2.45$\times 10^{-8}$
Table 3: Performance of various Riemannian methods on problem II (Model_Ob). $(n,k)$ | (40,8) | (50,10)
---|---|---
| Success | Time [s] | Iter. | Error | Success | Time [s] | Iter. | Error
RALM | 1 | 2.510 | 51 | 5.04$\times 10^{-7}$ | 0.95 | 4.727 | 64 | 4.94$\times 10^{-7}$
REPM${}_{\text{lqh}}$ | 0 | - | - | - | 0 | - | - | -
REPM${}_{\text{lse}}$ | 0 | - | - | - | 0 | - | - | -
RSQP | 0.65 | 8.618 | 5 | 2.30$\times 10^{-10}$ | 0.7 | 2.782$\times 10$ | 6 | 1.12$\times 10^{-10}$
RIPM | 1 | 3.791 | 22 | 5.62$\times 10^{-9}$ | 1 | 5.880 | 23 | 7.93$\times 10^{-9}$
$(n,k)$ | (60,12) | (70,14)
| Success | Time [s] | Iter. | Error | Success | Time [s] | Iter. | Error
RALM | 0.6 | 5.725 | 49 | 3.82$\times 10^{-7}$ | 0.6 | 8.223 | 52 | 3.85$\times 10^{-7}$
REPM${}_{\text{lqh}}$ | 0 | - | - | - | 0 | - | - | -
REPM${}_{\text{lse}}$ | 0 | - | - | - | 0 | - | - | -
RSQP | 0.7 | 4.446$\times 10$ | 5 | 1.17$\times 10^{-9}$ | 0.5 | 9.138$\times 10$ | 5 | 1.82$\times 10^{-9}$
RIPM | 1 | 7.134 | 23 | 9.69$\times 10^{-9}$ | 1 | 9.268 | 24 | 1.06$\times 10^{-8}$

##  9 Conclusions

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 first 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 efficiency of our
method.

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) objective $f\colon D\rightarrow\mathbb{R}$
defined on a geodesically convex subset $D\subset\mathcal{M}$. In contrast, in
(RCOP), we do not require any convexity. In practice, many convex functions
(in the Euclidean sense) are not geodesically convex on some interested
manifolds. For example, for any geodesically convex function defined on a
connected, compact Riemannian manifold (e.g., Stiefel manifold), it must be
constant [5, Corollary 11.10], which is not of interest in the field 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. 1.

RIPM generalizes EIPM from Euclidean space to general Riemannian manifolds.
EIPM is a special case of RIPM when $\mathcal{M}=\mathbb{R}^{n}$ or
$\mathbb{R}^{m\times n}$ in (RCOP).

  2. 2.

RIPM inherits all the advantages of Riemannian optimization. For example, we
can exploit the geometric structure of $\mathcal{M}$, which is usually
regarded as a set of constraints from the Euclidean viewpoint.

  3. 3.

Note that in both RIPM and EIPM, we have to solve the condensed Newton
equation (19) at each iteration. However, if the equality constraints can be
considered to be a manifold, RIPM can solve (19) with a _smaller_ order on
$T_{x}\mathcal{M}\times\mathbb{R}^{l}$. For example, the problem II (Model_St)
can be rewritten as

$\min_{X\in\mathbb{R}^{n\times k}}-2\operatorname{trace}(X^{T}C)\quad\text{
s.t. }X^{T}X=I_{k},\;X\geq 0.$

Here, Stiefel manifold is replaced by the equality constraints, i.e., we
define $h\colon\mathbb{R}^{n\times k}\to\operatorname{Sym}(k):X\mapsto
h(X):=X^{T}X-I_{k}$; and $\mathcal{M}=\mathbb{R}^{n\times k}$,
$l=\operatorname{dim}\operatorname{Sym}(k)=k(k+1)/2$ in (RCOP). Then, when we
apply EIPM, it requires us to solve (19) of order $nk+k(k+1)/2$. On the other
hand, if we apply RIPM to (Model_St), then (19) reduces to (21) since there
are only inequality constraints on $\mathcal{M}=\mathrm{St}(n,k)$. In this
case, we solve the equation of order $nk-k(k+1)/2$, i.e., the dimension of
$\mathrm{St}(n,k)$. Compared to EIPM, using RIPM reduces our dimensionality by
$k(k+1)$.

  4. 4.

RIPM can solve some problems that EIPM cannot. For example, the problem I
(NLRM) can be rewritten as

$\min_{X\in\mathbb{R}^{m\times n}}\|A-X\|_{\mathrm{F}}^{2}\quad\text{ s.t.
}\operatorname{rank}(X)=r,\;X\geq 0.$

Since the rank function, $X\mapsto\operatorname{rank}(X)$, is not even
continuous, we cannot apply EIPM.

##### Future Work I: Preconditioner for linear operator equation.

With regard to the complementary condition, $S_{k}^{-1}Z_{k}$ values display a
huge difference in magnitude as $k\to\infty$. The operator
$\Theta:=G_{x}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 to a problem that
does not have a matrix form. A possible alternative is to find a nonsingular
$\mathcal{P}$ such that the condition number of $\mathcal{P}^{-1}\mathcal{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 first
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 is done by monitoring only the KKT residuals; moreover, the
boundedness assumption (C2) of $\\{z_{k}\\}$ is too strong to hold in some
simple cases (see Wächter-Biegler effect [32]). It remains an important issue
to adapt more modern interior point methods to manifolds, although we may
encounter various difficulties 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.

Data availability The datasets generated during and/or analyzed during the
current study are available 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 qualifications 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, Birkhä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-Moré condition, in High-Performance Scientific 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 and applications, 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.
  * [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 with constraints, 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. Wä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, Pacific 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.