# Convergence and complexity of block majorization-minimization
for constrained block-Riemannian optimization

Yuchen Li  Department of Mathematics, University of Wisconsin-Madison, 480
Lincoln Dr., Madison, WI 53706, USA<EMAIL_ADDRESS>,  Laura Balzano
Department of Electrical Engineering and Computer Science, University of
Michigan, Ann Arbor, MI 48109, USA<EMAIL_ADDRESS>,  Deanna Needell
Department of Mathematics
University of California, Los Angeles, CA 90025, USA<EMAIL_ADDRESS>and
Hanbaek Lyu  Department of Mathematics, University of Wisconsin-Madison, 480
Lincoln Dr., Madison, WI 53706, USA<EMAIL_ADDRESS>
###### Abstract.

Block majorization-minimization (BMM) is a simple iterative algorithm for
nonconvex optimization that sequentially minimizes a majorizing surrogate of
the objective function in each block coordinate while the other block
coordinates are held fixed. We consider a family of BMM algorithms for
minimizing smooth nonconvex objectives, where each parameter block is
constrained within a subset of a Riemannian manifold. We establish that this
algorithm converges asymptotically to the set of stationary points, and
attains an $\epsilon$-stationary point within
$\widetilde{O}(\varepsilon^{-2})$ iterations. In particular, the assumptions
for our complexity results are completely Euclidean when the underlying
manifold is a product of Euclidean or Stiefel manifolds, although our analysis
makes explicit use of the Riemannian geometry. Our general analysis applies to
a wide range of algorithms with Riemannian constraints: Riemannian MM, block
projected gradient descent, optimistic likelihood estimation, geodesically
constrained subspace tracking, robust PCA, and Riemannian CP-dictionary-
learning. We experimentally validate that our algorithm converges faster than
standard Euclidean algorithms applied to the Riemannian setting.

###### Contents

  1. 1 Introduction
    1. 1.1 Related work
    2. 1.2 Our contribution
    3. 1.3 Organization
  2. 2 Preliminaries and Algorithm
    1. 2.1 Preliminaries on Riemannian geometry
    2. 2.2 Notations for block Riemannian optimization
    3. 2.3 Statement of algorithm
  3. 3 Statement of results
    1. 3.1 Optimality and complexity measures
    2. 3.2 Statement of results
  4. 4 Applications
    1. 4.1 Examples of manifolds
    2. 4.2 Euclidean block MM
    3. 4.3 Block Riemannian proximal updates on Hadamard manifolds
    4. 4.4 Block Riemannian/Euclidean proximal and prox-linear updates on Stiefel manifolds
    5. 4.5 Iteration complexity of MM on Stiefel manifolds in [16]
  5. 5 Stylized Applications
    1. 5.1 Geodesically constrained subspace tracking [18]
    2. 5.2 Optimistic likelihood under Fisher-Rao distance [49]
    3. 5.3 Riemannian CP-dictionary-learning
    4. 5.4 Robust PCA
  6. 6 Convergence analysis
    1. 6.1 Preliminary analysis
    2. 6.2 RBMM with Riemannian proximal surrogates
    3. 6.3 RBMM with Euclidean proximal surrogates
    4. 6.4 RBMM with $g$-smooth surrogates
    5. 6.5 Asymptotic stationarity
    6. 6.6 Rate of convergence
  7. A Background on Riemannian Optimization
  8. B Auxiliary lemmas
  9. C Details of Section 5
    1. C.1 Riemannian Hessian on fixed-rank manifold
    2. C.2 Details of Section 5.1
    3. C.3 Details of Section 5.2
    4. C.4 Details of Section 5.4

##  1\. Introduction

Optimization over Riemannian manifolds has a wide array of applications
ranging from the computation of linear algebraic quantities and factorizations
and problems with nonlinear differentiable constraints to the analysis of
shape space and automated learning [54, 35]. These applications arise either
because of implicit constraints on the problem to be optimized or because the
domain is naturally defined as a manifold. Although not nearly as abundant as
in the Euclidean optimization setting, there are various methods employed for
such optimization [11, 60, 10, 27]. Most of these approaches follow standard
optimization techniques by iteratively computing a descent direction and
taking a step in that direction along a geodesic. This computation is often
challenging in practice so other approaches alleviate this burden by utilizing
approximations or imposing additional assumptions [7, 17]. Here, we consider a
block majorization-minimization approach, based on the well-known
majorization-minimization (MM) method [38]. The MM method, as its name
suggests, consists of two main steps. The majorization step identifies an
upper bounding surrogate for the objective function and the minimization step
then minimizes that surrogate. The approach leads to iterations that decrease
the objective function, with the benefit that each minimization of a well-
chosen surrogate is computationally more efficient than the original problem.
Such an approach has also been applied in the context of unconstrained
manifold optimization with applications to robust sparse PCA on the Stiefel
manifold [16]. Here unconstrained manifold optimization means the constraint
set is only the manifold itself. When the underlying manifold is a Hadamard
manifold (i.e., complete, simply connected manifold with nonpositive sectional
curvature everywhere), then one can consider majorizing the objective function
by adding the Riemannian distance squared function. The resulting MM algorithm
corresponds to the Riemannian proximal point method [40, 13].

In this paper, we consider the minimization of a sufficiently smooth function
$f:\boldsymbol{\Theta}=\Theta^{(1)}\times\dots\times\Theta^{(m)}\rightarrow\mathbb{R}$.
Each constraint set $\Theta^{(i)}$ is a closed subset of a Riemannian manifold
$\mathcal{M}^{(i)}$. More precisely, we seek to solve

$\displaystyle\min_{\begin{subarray}{c}\boldsymbol{\theta}=[\theta^{(1)},\dots,\theta^{(m)}]\\\
\theta^{(i)}\in\Theta^{(i)}\subseteq\mathcal{M}^{(i)}\,\,\textup{for
$i=1,\dots,m$}\end{subarray}}f(\boldsymbol{\theta}).$ (1)

As the problem (1) is typically nonconvex, it is not always reasonable to
expect that an algorithm would converge to a globally optimal solution
starting from an arbitrary initialization. Instead, we aim to provide global
convergence (from arbitrary initialization) to stationary points. In some
problem classes, stationary points could be as good as global optimizers
practically as well as theoretically (see [44, 57]). Furthermore, obtaining
the iteration complexity of such algorithms is of importance both for
theoretical and practical purposes. In particular, one aims to bound the
worst-case number of iterations to achieve an $\varepsilon$-approximate
stationary point (defined appropriately, see Sec. 3.1).

In order to obtain a first-order optimal solution to (1), we consider various
Riemannian generalizations of the Block Majorization-Minimization (BMM)
algorithm in Euclidean space [33]. The high-level idea of BMM is that, in
order to minimize a multi-block objective, one can minimize a majorizing
surrogate of the objective in each block in a cyclic order. The algorithm we
study in the present work, which we call Riemannian Block Majorization-
Minimization (RBMM), can be stated in a high-level as follows (see Algorithm 1
for the full statement):

$\displaystyle\textbf{{RBMM}:}\,\begin{cases}&\textup{For $i=1,2,\dots,m$:}\\\
&g_{n}^{(i)}\leftarrow\left[\textup{Majorizing surrogate of $\theta\mapsto
f_{n}^{(i)}(\theta):=f\left(\theta_{n}^{(1)},\cdots,\theta_{n}^{(i-1)},\theta,\theta_{n-1}^{(i+1)},\cdots,\theta_{n-1}^{(m)}\right)$}\right]\\\
&\theta_{n}^{(i)}\in\operatorname*{arg\,min}_{\theta\in\Theta^{(i)}}g_{n}^{(i)}(\theta).\end{cases}$
(2)

In this work, we carefully analyze RBMM (2) in various settings and obtain
first-order optimality guarantees and iteration complexity. Moreover, the
connection between our RBMM and other existing algorithms is studied,
providing complexity results to some of them for the first time in the
literature.

As a preview, we provide a special case of Corollary 3.7 of our main results
concerning iteration complexity of RBMM on Stiefel manifolds (manifolds of
orthonormal frames, see Ex. 4.1). Note the term $L$-smooth in the following
corollary indicates the gradient of the function is $L$-Lipschitz continuous
(see Sec. 3.2 for comparison with geodesic smoothness). The
$\widetilde{O}(\cdot)$ notation is the variant of “big-O" notation that
ignores the logarithmic factors.

######  Corollary 1.1 (Complexity of RBMM on Stiefel manifolds).

Suppose we are minimizing an $L$-smooth function $f$ on the product of the
Stiefel manifolds using RBMM (2) with $L^{\prime}$-smooth surrogates with
quadratic majorization gap

$\displaystyle g_{n}^{(i)}(\theta)-f_{n}^{(i)}(\theta)\geq
c\lVert\theta-\theta_{n-1}^{(i)}\rVert^{2}$ (3)

for all $n\geq 1$ and $i=1,\dots,m$ for some $c>0$. Then the iterates
asymptotically converge to the set of stationary points and the algorithm has
iteration complexity of $\widetilde{O}(\varepsilon^{-2})$.

For instance, Euclidean block-proximal updates or Euclidean prox-linear
updates on the product of Stiefel manifolds have iteration complexity
$\widetilde{O}(\varepsilon^{-2})$. To the best of our knowledge, these type of
iteration complexity results for block Riemannian optimization methods are new
to the literature. Note that the conditions we need to check for applying
Corollary 1.1 (and its generalization Cor. 3.7) are completely Euclidean.
However, our proof of these corollaries incorporates Riemannian geometry in a
substantial manner.

###  1.1. Related work

Block majorization-minimization methods in the Euclidean setting for convex
problems are studied in [53] and [33] with general surrogates, in [59] with
prox-linear surrogates. On the other hand, the majorization-minimization
method on Riemannian manifolds is studied in [40] and [13] with proximal
surrogates on Hadamard manifolds, in [24] with tangent prox-linear surrogates
on Stiefel manifolds. These works under the Riemannian setting, however, only
consider the MM method in the one-block case.

Recently, a tangential type of block coordinate descent has been discussed in
[28]. There, the authors established a sublinear convergence rate for an
interesting block-wise Riemannian gradient descent. However, the retraction
considered there is restricted to the exponential map, which excludes many
commonly used retractions in the literature.

Most recently, block coordinate descent on smooth manifolds has been studied
in [51]. The authors studied three types of block coordinate methods on
manifolds, namely block majorization-minimization (BMM), block coordinate
descent (BCD), and block Riemannian gradient descent (BRGD), as well as the
blending of the latter two. Asymptotic convergence to stationary points of BCD
and BMM is guaranteed by Theorem 1 and Theorem 2, respectively. Furthermore,
the convergence rate of BRGD, as well as its blending with BCD, is shown in
their Theorem 3 and Theorem 4, where the manifold is assumed to be a compact
smooth submanifold of $\mathbb{R}^{n}$. It is important to note, however, that
convergence (rate) results are applicable to unconstrained problems with exact
computations, i.e. an exact minimizer is needed for each sub-problem.
Additionally, the rate of convergence of BMM remains unaddressed.

In Table 1, we provide a summary of the aforementioned related work, along
with the details.

Methods | Manifold | Objective | Constraints | Blocks | Complexity | $\begin{matrix}\text{Inexact}\\\ \text{comp.}\end{matrix}$
---|---|---|---|---|---|---
Euclidean BMM [53] |  |  |  |  |  |
Block Prox-linear [59] | Euclidean | convex | convex | many | $\widetilde{O}(\varepsilon^{-1})$ | ✗
Block PGD [19] |  |  |  |  |  |
Riemannian prox. [40] | Hadamard |  $g$-convex |  $g$-convex | 1 | - | ✗
Riemannian prox. [13] | Hadamard |  $g$-convex |  $g$-convex | 1 | $\widetilde{O}(\varepsilon^{-1})$ | ✗
Riemannian Prox-linear | Riemannian & | non-convex & | N/A | 1 | $\widetilde{O}(\varepsilon^{-2})$ | ✗
(linear search)[24] | Compact |  smooth†
Block Riemannian GD | Riemannian | non-convex | N/A | many | $\widetilde{O}(\varepsilon^{-2})$ | ✗
(exponential map)[28]
BMM on manifolds[51] | Riemannian | non-convex | N/A | many | - | ✗
RBMM (Ours) with surrogates: |  |  |  |  |  |
$g$-smooth (Thm. 3.6) | Riemannian | non-convex |  $g$-convex | many | $\widetilde{O}(\varepsilon^{-4})^{\dagger}$ | ✓
Riemannian proximal | Riemannian | non-convex | $g$-convex | many | $\widetilde{O}(\varepsilon^{-2})$ | ✓
(Thm. 3.5)
Euclidean proximal | Riemannian | non-convex |  $g$-convex & | many | $\widetilde{O}(\varepsilon^{-2})$ | ✓
(Thm. 3.5) |  $\subseteq$ Euclidean | compact
Smooth (Cor. 3.7) | Euclidean/ | non-convex & | convex/ | many | $\widetilde{O}(\varepsilon^{-2})$ | ✓
Stiefel |  smooth† |  $g$-convex
Table 1. Our main contributions and comparison to existing results.
“$g$-smooth" means being geodesically smooth with respect to the geometry of
the underlying manifold, and “smooth” means being smooth with respect to the
Euclidean geometry. “$g$-convex” means geodesic convexity of subsets of
manifolds. $\widetilde{O}(\cdot)$ notations means big-$O$ up to logarithmic
factors. The objective function marked by “smooth†” only needs to be smooth in
the Euclidean sense; In all other cases, it is required to be $g$-smooth with
respect to the underlying manifold. The last column shows whether the method
allows the inexact solution to a subproblem, i.e. the robustness under inexact
computation. Details of comparison to known results can be found in Section 4.
The complexity of $\widetilde{O}(\varepsilon^{-4})$ in the first row of RBMM
can be improved to $\widetilde{O}(\varepsilon^{-2})$ under some additional
assumption, see Theorem 3.6.

###  1.2. Our contribution

In this work, we thoroughly analyze RBMM (2) and obtain asymptotic convergence
to stationary points and iteration complexity. The novelty of this work,
compared to the aforementioned related work, lies especially in the following
three aspects:

(1) (Constrained optimization) RBMM is applicable to constrained optimization
problems on manifolds. Here, constrained optimization on manifolds means we
allow the constraint set $\boldsymbol{\Theta}$ of the optimization problem to
be a closed subset of the manifold, i.e.
$\boldsymbol{\Theta}\subseteq\mathcal{M}$, which is not necessarily the entire
manifold.

(2) (Iteration complexity) While asymptotic convergence of certain RBMM
variants is established in existing literature, the rate of convergences
remains unknown. Our work fills this gap by deriving the iteration complexity
of RBMM. See Theorems 3.5 and 3.6.

(3) (Robustness) RBMM is robust in the face of inexact computation of
optimization sub-problems. See (A0)(ii).

RBMM, as a general Riemannian block optimization framework, entails many
classical algorithms including Euclidean block MM, proximal updates on
Hadamard manifolds, and MM methods on Stiefel manifolds. We apply our results
to various stylized applications such as geodesically constrained subspace
tracking, optimistic likelihood under Fisher-Rao distance, Riemannian CP-
dictionary-learning, and robust PCA and obtain the following results:

(4) Asymptotic convergence and complexity of $\widetilde{O}(\varepsilon^{-2})$
for Euclidean block proximal methods when the manifolds are embedded in
Euclidean spaces and the constraint sets on the manifolds are compact and
$g$-convex. See Theorem 3.5.

(5) Asymptotic convergence and complexity of $\widetilde{O}(\varepsilon^{-2})$
for block MM with Euclidean $L$-smooth surrogate on the product of Stiefel and
Euclidean manifolds with $g$-convex constraints. See Corollary 3.7.

(6) A Euclidean-regularized version of MM with linear surrogates on the
Stiefel manifold is proposed, with guarantees on asymptotic convergence and
complexity of $\widetilde{O}(\varepsilon^{-2})$, which can be applied to
geodesically constrained subspace tracking problems [18]. See Corollary 4.11
and Corollary 5.1.

Table 1 presents a summary of related work, along with our main theoretical
contributions.

###  1.3. Organization

The paper is organized as follows. We introduce preliminaries and notations in
Section 2.1 and 2.2. In Section 2.3, we give a precise statement of the RBMM
algorithm. We detail the standing assumptions and state our main results in
Section 3. In Section 4, we present some applications and classical algorithms
as special cases of our general results. We present some stylized applications
of our results in Section 5. We prove the convergence results of RBMM,
Algorithm 1, throughout Section 6. Figure 1 provides a structure diagram of
the present paper.

Figure 1.  Structure of the present paper.

##  2\. Preliminaries and Algorithm

###  2.1. Preliminaries on Riemannian geometry

The notations throughout this paper are consistent with the common literature,
see e.g. [7] and [17]. For background knowledge on Riemannian geometry, we
refer the readers to [56], [37], [25], and [32]. In this section, we provide a
brief introduction to the notations used in our paper, with further details
provided in the Appendix A.

A Riemannian manifold $\mathcal{M}$ is a manifold endowed with a Riemannian
metric
$\left(\eta,\xi\right)\mapsto\left\langle\eta,\xi\right\rangle_{x}\in\mathbb{R}$,
where $\eta$ and $\xi$ are tangent vectors in the tangent space
$T_{x}\mathcal{M}$ (also denoted as $T_{x}$ when it is clear from context) of
$\mathcal{M}$ at $x$. This inner product in the tangent space is also denoted
by $\langle\cdot,\cdot\rangle$ for convenience when the subscript is clear
from the context. The induced norm on the tangent space is denoted by
$\|\cdot\|_{x}$ or $\|\cdot\|$. The Riemannian gradient of a smooth function
$f:\mathcal{M}\rightarrow\mathbb{R}$ at $x$ is defined as the unique tangent
vector, $\operatorname{grad}f(x)\in T_{x}\mathcal{M}$, such that
$\langle\operatorname{grad}f(x),\xi_{x}\rangle=\mathrm{D}f(x)\left[\xi_{x}\right],\forall\xi_{x}\in
T_{x}\mathcal{M}$, where $\mathrm{D}f(x)\left[\xi_{x}\right]$ is the
differential of $f$ at point $x$ along the direction $\xi_{x}$. The geodesic
distance between $x,y\in\mathcal{M}$ is denoted by $d_{\mathcal{M}}(x,y)$ or
$d(x,y)$ when it is clear from the context.

A retraction on a manifold $\mathcal{M}$ is a locally defined smooth mapping
$\operatorname{\textup{Rtr}}$ from the tangent bundle $T\mathcal{M}$ to
$\mathcal{M}$ with the following properties.

(i)

For each $x\in\mathcal{M}$, let $\operatorname{r_{\textup{Rtr}}}(x)>0$ be the
‘retraction radius’ such that the restriction
$\operatorname{\textup{Rtr}}_{x}:T_{x}\mathcal{M}\rightarrow\mathcal{M}$ of
$\operatorname{\textup{Rtr}}$ to $T_{x}\mathcal{M}$ is well-defined in a ball
of radius $\operatorname{r_{\textup{Rtr}}}(x)$ around the origin
$\mathbf{0}=\mathbf{0}_{x}$.

(ii)

$\operatorname{\textup{Rtr}}_{x}(\mathbf{0})=x$; The differential of
$\operatorname{\textup{Rtr}}_{x}$ at $\mathbf{0}$,
$D\operatorname{\textup{Rtr}}_{x}(\mathbf{0})$, is the identity map on
$T_{x}\mathcal{M}$.

For each $x\in\mathcal{M}$ and $\eta\in T_{x}\mathcal{M}$, the retraction
curve $t\mapsto\operatorname{\textup{Rtr}}_{x}(t\eta)$ agrees up to first
order with geodesics passing through $x$ with velocity $\eta$ around $x$.
Retractions provide a way to lift a function
$g:\mathcal{M}\rightarrow\mathbb{R}$ onto the tangent space $T_{x}\mathcal{M}$
via its pullback
$\hat{g}:=g\circ\operatorname{\textup{Rtr}}_{x}:T_{x}\rightarrow\mathbb{R}$.
We use this construction to lift an upper-bounding surrogate defined on the
manifold onto the tangent spaces in Algorithm 1. Note that for all $\eta\in
T_{x}\mathcal{M}$,

$\displaystyle\langle\nabla\hat{g}(\mathbf{0}),\,\eta\rangle=D\hat{g}(\mathbf{0})[\eta]=Dg(x)[D\operatorname{\textup{Rtr}}_{x}(\mathbf{0})[\eta]]=Dg(x)[\eta]=\langle\operatorname{grad}g(x),\,\eta\rangle.$
(4)

If for all $x\in\mathcal{M}$ and $\eta\in T_{x}\mathcal{M}$, the retraction
curve $t\mapsto\operatorname{\textup{Rtr}}_{t\eta}$ coincides with the
geodesic curve $t\mapsto\gamma(t)$ with $\gamma(0)=x$ and
$\gamma^{\prime}(0)=\eta$ whenever $\lVert
t\eta\rVert\leq\operatorname{r_{\textup{Exp}}}(x)$ for some constant
$\operatorname{r_{\textup{Exp}}}(x)>0$, then the retraction
$\operatorname{\textup{Rtr}}$ is called the exponential map and denoted as
$\operatorname{Exp}$. Note that the exponential map is defined as the solution
of a nonlinear ordinary differential equation. While every Riemannian manifold
admits the exponential map, its computation is often challenging. Retractions
provide computationally efficient alternatives to exponential maps. Some
typical choices of retractions are
$\operatorname{\textup{Rtr}}_{x}(\eta)=x+\eta$ on Euclidean spaces and
$\operatorname{\textup{Rtr}}_{x}(\eta)=\frac{x+\eta}{\lVert x+\eta\rVert}$ on
spheres. See [7, Sec. 4.1] for more examples of retractions. If the
exponential map is defined on the entire tangent bundle (i.e.,
$\operatorname{r_{\textup{Exp}}}(x)=\infty$ for all $x\in\mathcal{M}$), then
we say $\mathcal{M}$ is (geodesically) complete. Note by definition of
exponential map we have $\|\eta\|=d(x,y)$ whenever
$d(x,y)\leq\operatorname{r_{\textup{Exp}}}(x)$ and
$\operatorname{Exp}_{x}(\eta)=y$.

A Riemannian manifold is locally diffeomorphic to its tangent spaces, so it
resembles the Euclidean space within a small metric ball around each point.
Accordingly, there are several notions of radius functions
$r:\mathcal{M}\rightarrow[0,\infty)$, including the injectivity and the
convexity radii. For $x\in\mathcal{M}$, consider the open ball
$B(x,r)=\\{\eta\in T_{x}\mathcal{M}:\langle\eta,\eta\rangle<r\\}\subseteq
T_{x}\mathcal{M}$; the injectivity radius of $\mathcal{M}$ at $x$, denoted as
$\operatorname{r_{\textup{inj}}}(x)$, is the supremum of values of $r$ such
that $\operatorname{Exp}_{x}$ defines a diffeomorphism from $B(x,r)$ to its
image on $\mathcal{M}$. Thus, we can also define the inverse exponential map
from $\mathcal{M}$ to $T_{x}\mathcal{M}$, denoted by
$\operatorname{Exp}^{-1}_{x}(\cdot)$, within the injectivity radius. A set
$C\subseteq\mathcal{M}$ is called (geodesically) strongly convex if for any
$x$ and $y$ in $C$, there is a unique minimal geodesic $\gamma$ in
$\mathcal{M}$ joining $x$ and $y$, and $\gamma$ is contained in $C$. In
particular, compact manifolds have uniformly positive injectivity radius (see
[22, Thm. III.2.3]). Furthermore, it is worth noting that Hadamard manifolds
(complete and simply connected manifolds with non-positive curvature), which
include Euclidean space, hyperbolic space, and manifolds of positive definite
matrices, also have uniformly positive injectivity radius ([4], [56, Theorem
4.1, p.221]). When the injectivity radius is uniformly positive, there exists
a retraction with a uniformly positive retraction radius (e.g., the
exponential map).

For a subset $\boldsymbol{\Theta}\subseteq\mathcal{M}$ and
$x\in\boldsymbol{\Theta}$, define the lifted constraint set
$T_{x}^{*}\mathcal{M}$ as

$\displaystyle T^{*}_{x}\mathcal{M}$ $\displaystyle:=\\{u\in
T_{x}\mathcal{M}\,|\,\textup{$\operatorname{\textup{Rtr}}_{x}(u)=x^{\prime}$
for some $x^{\prime}\in\boldsymbol{\Theta}$ with $d(x,x^{\prime})\leq
r_{0}/2$}\\},$ (5)

where $r_{0}$ is the lower bound of the injectivity radius (see (A0) in
Section 3). When we use exponential map $\operatorname{Exp}$ as the retraction
in (5), one can think of the set $T^{*}_{x}\mathcal{M}$ as the ‘lift’ of the
constraint set $\boldsymbol{\Theta}$ onto the tangent space $T_{x}\mathcal{M}$
in the sense that if $\boldsymbol{\Theta}$ is contained in the metric ball of
radius $\operatorname{r_{\textup{inj}}}(x)$ centered at $x$, then
$T^{*}_{x}\mathcal{M}$ equals the inverse image
$\operatorname{Exp}^{-1}_{x}(\boldsymbol{\Theta})$ of $\boldsymbol{\Theta}$
under $\operatorname{Exp}_{x}$. In particular, this formula holds for all
subset $\boldsymbol{\Theta}$ of $\mathcal{M}$ if $\mathcal{M}$ is a complete
Riemannian manifold since $\operatorname{Exp}_{x}$ is defined on the entire
$T_{x}\mathcal{M}$. If $\mathcal{M}$ is a Euclidean space and
$\boldsymbol{\Theta}$ is a convex subset of it, then
$T^{*}_{x}\mathcal{M}=\boldsymbol{\Theta}$. Lastly, we note that when
$\boldsymbol{\Theta}$ is strongly convex in $\mathcal{M}$, the set
$T_{x}^{*}\mathcal{M}$ above is locally defined near $x$. That is, we can
replace the injectivity radius $r_{0}$ in (5) by any constant
$\delta\in(0,r_{0})$.

On Riemannian manifolds, parallel transport provides a way to transport a
vector along a smooth curve. For each smooth curve
$\gamma:[0,1]\rightarrow\mathcal{M}$, denote the parallel transport along
$\gamma$ from $x=\gamma(0)$ to $y=\gamma(t)$ by $\Gamma^{\gamma}_{x\rightarrow
y}$. If $\gamma$ is clear from the context (e.g., the unique distance-
minimizing geodesic from $x$ to $y$), then we also write $\Gamma_{x\rightarrow
y}^{\gamma}=\Gamma_{x\rightarrow y}$. Intuitively, a tangent vector $\eta\in
T_{x}\mathcal{M}$ at $x$ of $\gamma$ is still a tangent vector
$\Gamma(\gamma)\eta\in T_{y}\mathcal{M}$ of $\gamma$ at $y$. Recall that, one
important property of parallel transport on a Riemannian manifold is that it
preserves the inner product, i.e.

$\displaystyle\left\langle\Gamma^{\gamma}_{\gamma(0)\rightarrow\gamma(t)}\,(\xi),\,\Gamma^{\gamma}_{\gamma(0)\rightarrow\gamma(t)}\,(\zeta)\right\rangle_{\gamma(t)}=\left\langle\xi,\,\zeta\right\rangle_{\gamma(0)}\quad\textup{for
all $t\in[0,1]$, $\xi,\zeta\in T_{\gamma(0)}$}.$ (6)

See Figure 2 for an illustration.

Figure 2. Illustration of (left) retraction and (right) parallel transport.

###  2.2. Notations for block Riemannian optimization

In (1), we are interested in minimizing an objective function $f$ within the
product parameter space
$\boldsymbol{\Theta}=\Theta^{(1)}\times\dots\times\Theta^{(m)}$, where each
constraint set $\Theta^{(i)}$ is a subset of a Riemannian manifold
$\mathcal{M}^{(i)}$. It will be convenient to introduce the following
notations: For $\boldsymbol{\theta}=[\theta^{(1)},\dots,\theta^{(m)}]$,

$\displaystyle\operatorname{grad}_{i}f(\boldsymbol{\theta})$
$\displaystyle:=\textup{Riemmanian gradient of $\theta\mapsto
f(\theta^{(1)},\dots,\theta^{(i-1)},\theta,\theta^{(i+1)},\dots,\theta^{(m)})$},$
(7) $\displaystyle\operatorname{grad}f(\boldsymbol{\theta})$
$\displaystyle:=[\operatorname{grad}_{1}f(\boldsymbol{\theta}),\dots,\operatorname{grad}_{m}f(\boldsymbol{\theta})],$
(8) $\displaystyle d(\mathbf{x},\mathbf{y})$
$\displaystyle:=\sqrt{\sum_{i=1}^{m}d(x^{(i)},y^{(i)})^{2}}\quad\textup{for
$\mathbf{x}=(x^{(1)},\dots,x^{(m)}),\,\mathbf{y}=(y^{(1)},\dots,y^{(m)})\in\prod_{i=1}^{m}\mathcal{M}^{(i)}$}.$
(9)

Note that if we endow the product $\prod_{i=1}^{m}\mathcal{M}^{(i)}$ of the
manifolds a joint Riemannian structure, then we can interpret
$\operatorname{grad}f(\boldsymbol{\theta})$ above as the Riemannian gradient
at $\boldsymbol{\theta}$ with respect to that joint Riemannian structure.
However, we do not explicitly introduce or use such a product manifold
structure in the manuscript.

Throughout this paper, we let $(\boldsymbol{\theta}_{n})_{n\geq 1}$ denote an
output of Algorithm 1 and write
$\boldsymbol{\theta}_{n}=[\theta_{n}^{(1)},\dots,\theta_{n}^{(m)}]$ for each
$n\geq 1$. For each $n\geq 1$ and $i=1,\dots,m$, denote

$\displaystyle f_{n}^{(i)}:\theta\mapsto
f(\theta_{n}^{(1)},\dots,\theta_{n}^{(i-1)},\theta,\theta_{n-1}^{(i+1)},\dots,\theta_{n-1}^{(m)}),$
(10)

which we will refer to as the $i$th marginal objective function at iteration
$n$.

###  2.3. Statement of algorithm

Below in Algorithm 1, we give a precise statement of the RBMM algorithm we
stated in high-level at (2). We first define majorizing surrogate functions
defined on Riemannian manifolds.

######  Definition 2.1 (Majorizing surrogates on Riemannian manifolds).

Fix a function $h:\mathcal{M}\rightarrow\mathbb{R}$ and a point
$\boldsymbol{\theta}\in\mathcal{M}$, where $\mathcal{M}$ is a Riemannian
manifold. A function $g:\mathcal{M}\rightarrow\mathbb{R}$ is a majorizing
surrogate of $h$ at $\boldsymbol{\theta}$ if

$\displaystyle g(x)\geq h(x)\quad\textup{for all
$x\in\mathcal{M}$}\quad\textup{and}\quad
g(\boldsymbol{\theta})=h(\boldsymbol{\theta}).$ (11)

As mentioned before, the high-level idea is the following. In order to update
the $i$th block of the parameter $\theta^{(i)}_{n}$ at iteration $n$, we use
the RMM we described in the introduction.

Algorithm 1 Riemannian Block Majorization-Minimization (RBMM)

1:Input:
$\boldsymbol{\theta}_{0}=(\theta_{0}^{(1)},\cdots,\theta_{0}^{(m)})\in\Theta^{(1)}\times\cdots\times\Theta^{(m)}$
(initial estimate); $N$ (number of iterations); $\alpha\in[0,1]$ (step size)

2: for $n=1,\dots,N$ do:

3: Update estimate
$\boldsymbol{\theta}_{n}=[\theta_{n}^{(1)},\cdots,\theta_{n}^{(m)}]$ by

4: For $i=1,\cdots,m$ do:

5: $\displaystyle
f_{n}^{(i)}(\cdot):=f\left(\theta_{n}^{(1)},\cdots,\theta_{n}^{(i-1)},\,\cdot\,\,,\theta_{n-1}^{(i+1)},\cdots,\theta_{n-1}^{(m)}\right):\mathcal{M}^{(i)}\rightarrow\mathbb{R}\qquad\textup{($\triangleright$
marginal objective function)}$

6:
$\displaystyle\qquad\begin{cases}&g_{n}^{(i)}\leftarrow\left[\textup{Majorizing
surrogate of $f_{n}^{(i)}$ at $\theta_{n-1}^{(i)}$}\right]\\\
&\theta_{n}^{(i)}\in\operatorname*{arg\,min}_{\theta\in\Theta^{(i)}}g_{n}^{(i)}(\theta)\end{cases}$
(12)

7: end for

8: end for

9:output: $\boldsymbol{\theta}_{N}$

In Algorithm 1, the majorizing surrogate $g_{n}^{(i)}$ at each iteration $n$
for each block $i$ is chosen so that

(1) (Majorization) $g_{n}^{(i)}(x)-f_{n}^{(i)}(x)\geq 0$ for all
$x\in\mathcal{M}^{(i)}$;

(2) (Sharpness)
$g_{n}^{(i)}(\theta_{n-1}^{(i)})=f_{n}^{(i)}(\theta_{n-1}^{(i)})$.

##  3\. Statement of results

In this section, we state our main results concerning the convergence and
complexity of our RBMM algorithm (Alg. 1) for the constrained block Riemannian
optimization problem in (1). Figure 1 provides a structure diagram of the main
results and assumptions.

###  3.1. Optimality and complexity measures

For iterative algorithms, first-order optimality conditions may hardly be
satisfied exactly in a finite number of iterations, so it is more important to
know how the worst-case number of iterations required to achieve an
$\varepsilon$-approximate solution scales with the desired precision
$\varepsilon$.

More precisely, for the multi-block problem (1), we say
$\boldsymbol{\theta}_{*}=[\theta_{*}^{(1)},\dots,\theta_{*}^{(m)}]\in\boldsymbol{\Theta}$
is an $\varepsilon$-stationary point of $f$ over $\boldsymbol{\Theta}$ if

$\displaystyle\sum_{i=1}^{m}\,\left(-\inf_{u\in
T_{\theta_{*}^{(i)}}^{*},\|u\|\leq
1}\left\langle\operatorname{grad}_{i}f(\boldsymbol{\theta}_{*}),\,\frac{u}{\hat{r}}\right\rangle\right)\leq\varepsilon,$
(13)

where $\hat{r}=\min\\{r_{0},1\\}$, i.e. the minimum of $1$ and the lower bound
of injectivity radius.

In the Euclidean setting, (13) reduces to the optimality measure used in [43]
for constrained nonconvex smooth optimization and is also equivalent to [48,
Def. 1] for smooth objectives.

In the unconstrained block-Riemannian setting where
$\Theta^{(i)}=\mathcal{M}^{(i)}$ for $i=1,\dots,m$, the above equation (13)
becomes

$\displaystyle\sum_{i=1}^{m}\,\lVert\operatorname{grad}_{i}f(\boldsymbol{\theta}_{*})\rVert\leq\varepsilon.$
(14)

In the case of single-block $m=1$, the above is the standard definition of
$\varepsilon$-stationary points for unconstrained Riemannian optimization
problems.

Next, for each $\varepsilon>0$ we define the (worst-case) iteration complexity
$N_{\varepsilon}$ of an algorithm computing $(\boldsymbol{\theta}_{n})_{n\geq
1}$ for solving (1) as

$\displaystyle
N_{\varepsilon}:=\sup_{\boldsymbol{\theta}_{0}\in\boldsymbol{\Theta}}\,\inf\,\\{n\geq
1\,|\,\text{$\boldsymbol{\theta}_{n}$ is an $\varepsilon$-approximate
stationary point of $f$ over $\boldsymbol{\Theta}$}\\},$ (15)

where $(\boldsymbol{\theta}_{n})_{n\geq 0}$ is a sequence of estimates
produced by the algorithm with an initial estimate $\boldsymbol{\theta}_{0}$.
Note that $N_{\varepsilon}$ gives the worst-case bound on the number of
iterations for an algorithm to achieve an $\varepsilon$-approximate solution
due to the supremum over the initialization $\boldsymbol{\theta}_{0}$ in (15).

###  3.2. Statement of results

Here we state our main convergence results of solving the minimization problem
(1) using Algorithm 1. First, we introduce a Riemannian counterpart of a
smoothness property of a function $F:\mathcal{M}\rightarrow\mathbb{R}$. In the
Euclidean setting, The function $F$ is $L$-smooth if its gradient $\nabla F$
is $L$-Lipschitz continuous. In the Riemannian case, Riemannian gradients
$\operatorname{grad}F(x)$ and $\operatorname{grad}F(y)$ at two base points
$x,y\in\mathcal{M}$ live in different tangent spaces $T_{x}$ and $T_{y}$, so
they have to be compared using a parallel transport (see Figure 2 for an
illustration). We extend this notion of smoothness to the setting where the
function $F$ is defined on the product of Riemannian manifolds.

######  Definition 3.1 (Geodesic smoothness).

A function $F:\prod_{i=1}^{m}\mathcal{M}^{(i)}\to\mathbb{R}$ is geodesically
smooth ($g$-smooth in short) with parameter $L>0$ if $F$ is block-wise
continuously differentiable and for each
$\mathbf{x}=(x^{(1)},\dots,x^{(m)}),\mathbf{y}=(y^{(1)},\dots,y^{(m)})\in\prod_{i=1}^{m}\mathcal{M}^{(i)}$
where there exists a minimizing geodesic joining $x^{(i)}$ and $y^{(i)}$ for
each $i=1,\dots,m$,

$\left\lVert\operatorname{grad}_{i}F(\mathbf{x})-\Gamma_{y^{(i)}\rightarrow
x^{(i)}}(\operatorname{grad}_{i}F(\mathbf{y}))\right\rVert\leq\frac{L}{m}d(\mathbf{x},\mathbf{y}),$
(16)

where $\Gamma_{y^{(i)}\rightarrow x^{(i)}}$ is the parallel transport along a
distance-minimizing geodesic joining $x^{(i)}$ and $y^{(i)}$ in
$\mathcal{M}^{(i)}$, and $d(\mathbf{x},\mathbf{y})$ is defined in (9).

An important consequence of the $g$-smoothness is the following quadratic
bound on first-order approximation:

$\left|F(y)-F(x)-\left\langle\operatorname{grad}F(x),\gamma^{\prime}(0)\right\rangle_{x}\right|\leq\frac{L}{2}d^{2}(x,y),$
(17)

where $\gamma$ is any distance-minimizing geodesic in $\mathcal{M}$ from $x$
to $y$ and $d(x,y)$ is the Riemannian distance between $x$ and $y$. See
Appendix B.1 for the proof. Throughout this paper, the term “$g$-smooth”
denotes geodesic smooth (Def 3.1), while “$L$-smooth” indicates the function
$L$-smooth in Euclidean sense, meaning the Euclidean gradient of the function
is $L$-Lipschitz continuous.

We start by stating some general assumptions. We allow inexact computation of
the solution to the minimization sub-problems in Algorithm 1. This is
practical since minimizing the surrogates on the manifold (possibly with
additional constraints) may not always be exactly solvable. To be precise, for
each $n\geq 1$, we define the optimality gap $\Delta_{n}$ by

$\Delta_{n}=\Delta_{n}\left(\boldsymbol{\theta}_{0}\right):=\max_{1\leq i\leq
m}\left(g_{n}^{(i)}(\theta_{n}^{(i)})-\inf_{\theta\in\boldsymbol{\Theta}^{(i)}}g_{n}^{(i)}(\theta)\right)$
(18)

For the convergence analysis to hold, we require that the optimal gaps decay
fast enough so that they are summable, stated as Assumption (A0).

######  (A0).

For RBMM (Alg. 1), we make the following assumptions:

(i)

(Objective) The objective
$f:\boldsymbol{\Theta}=\prod_{i=1}^{m}\Theta^{(i)}\rightarrow\mathbb{R}$ is
$g$-smooth with some parameter $L_{f}>0$ and its values are uniformly lower
bounded by some $f^{*}\in\mathbb{R}$. Furthermore, the sublevel sets
$f^{-1}((-\infty,a))=$
$\\{\boldsymbol{\theta}\in\boldsymbol{\Theta}:f(\boldsymbol{\theta})\leq a\\}$
are compact for each $a\in\mathbb{R}$.

(ii)

(Inexact computation) The optimality gaps $\Delta_{n}$ in (18) are summable,
that is, $\sum_{n=1}^{\infty}\Delta_{n}<\infty$. Furthermore, let
$\theta_{n}^{(i\star)}$ be an exact solution of the minimization step in
(LABEL:eq:RBMM_MmMm) and let $\theta_{n}^{(i)}$ be the inexact output. Then
for all $i=1,\dots,m$,

$\displaystyle d(\theta_{n}^{(i\star)},\theta_{n}^{(i)})=o(1).$ (19)

If the surrogate function $g^{(i)}_{n}$ is geodesically strongly convex (see
Definition 6.5), then (19) is a direct consequence of the summability of the
optimality gaps. In particular, this is the case for Riemannian proximal
surrogates on Hadamard manifolds as a special case of (A1)(i-rp) (see Prop.
6.7).

Next, we impose the following conditions on the underlying manifolds,
majorizing surrogates, and constraint sets.

######  (A1).

For the manifolds and surrogates, at least one of the following holds:

(i-gs)

($g$-smooth surrogates) Each surrogate
$g_{n}^{(i)}:\mathcal{M}^{(i)}\rightarrow\mathbb{R}$ is $g$-smooth with some
parameter $L_{g}\geq 0$ for all $n\geq 1$ and $i=1,\dots,m$.

(i-rp)

(Riemannian Proximal surrogates) Each surrogate $g_{n}^{(i)}$ is a Riemannian
proximal surrogate; that is, for each $n\geq 1$ and some constant
$\lambda_{n}\geq L_{f}$ with $\lambda_{n}=O(1)$,

$\displaystyle
g_{n}^{(i)}(\theta)=f_{n}^{(i)}(\theta)+\frac{\lambda_{n}}{2}d^{2}(\theta,\theta_{n-1}^{(i)}).$
(20)

(i-ep)

(Euclidean Proximal surrogates) Each $\mathcal{M}^{(i)}$ is an embedded
submanifold in a Euclidean space and the constraint set $\Theta^{(i)}$ is
compact. Each surrogate $g_{n}^{(i)}$ is a Euclidean proximal surrogate, that
is, for each $n\geq 1$ and some constant $\lambda_{n}\geq L_{f}$ with
$\lambda_{n}=O(1)$,

$\displaystyle
g_{n}^{(i)}(\theta)=f_{n}^{(i)}(\theta)+\frac{\lambda_{n}}{2}\lVert\theta-\theta_{n-1}^{(i)}\rVert^{2}.$
(21)

Furthermore, we require the following for the constraint sets:

(ii)

For each $i=1,\dots,m$, there exists a uniform lower bound $r_{0}>0$ for the
injectivity radius $\operatorname{r_{\textup{inj}}}(x)$ over
$x\in\Theta^{(i)}$, i.e. $\operatorname{r_{\textup{inj}}}(x)\geq r_{0}$ for
all $x\in\Theta^{(i)}$.

When the manifolds $\mathcal{M}^{(i)}$ are non-Euclidean, constructing
$g$-smooth surrogates may not always be easy. However, for Stiefel manifolds
(see Ex. 4.1), Euclidean smoothness and $g$-smoothness are essentially
equivalent (see Lem. 4.9) so standard smooth surrogates in the Euclidean space
(e.g, proximal (21) and prox-linear (36)) verifies (A1)(i-gs). This fact is
crucially used in the proof of Corollaries 1.1 and 3.7.

Note that both the Riemannian and the Euclidean proximal surrogates in
(20)-(21) are not necessarily $g$-smooth (see Section 6.2 for details). Our
analysis of RBMM with these proximal surrogates uses explicit forms of the
Riemannian/Euclidean gradient of the squared geodesic/Euclidean distance
function (see Prop. 6.8) instead of $g$-smoothness of the surrogates.

The surrogates in (A1)(i-rp) and (A1)(i-ep) are identical when the manifold
$\mathcal{M}^{(i)}$ is the Euclidean space. For non-Euclidean manifolds, the
difference between these two surrogates could be significant since the
Riemannian proximal surrogates in (A1)(i-rp) use the geometry of the non-
Euclidean manifold $\mathcal{M}^{(i)}$ while the Euclidean proximal surrogates
in (A1)(i-ep) use the geometry of the ambient Euclidean space. For example, in
the case of a sphere, (A1)(i-rp) uses the arc length of great circles while
(A1)(i-ep) uses length of line segments. While the Riemannian proximal
surrogates in (A1)(i-rp) can handle non-compact manifolds without compact
constraints (e.g., see Sec. 5.2), the Euclidean proximal surrogates in
(A1)(i-ep), if applicable, are computationally simpler to handle (e.g., ease
of differentiation). We also remark that (A1)(i-ep) holds for unconstrained
optimization problems on compact manifolds embedded in Euclidean space, i.e.
$\Theta^{(i)}=\mathcal{M}^{(i)}$ and $\mathcal{M}^{(i)}$ is a compact embedded
submanifold of a Euclidean space (e.g., Stiefel manifolds).

We now state our first main result, concerning the asymptotic convergence of
RBMM for the constrained single/two-block block Riemannian optimization
problem (1).

######  Theorem 3.2 (Asymptotic convergence to stationary points; one or two
blocks).

Let $f$ denote the objective function in (1) with $m=1$ or $2$. Let
$(\boldsymbol{\theta}_{n})_{n\geq 0}$ be an output of Algorithm 1 under (A0)
and (A1)(i-gs),(ii). Further assume that the constraint set $\Theta^{(i)}$ of
$\mathcal{M}^{(i)}$ is strongly $g$-convex for $i=1,\dots,m$. Then every limit
point of $(\boldsymbol{\theta}_{n})_{n\geq 0}$ is a stationary point of $f$
over $\boldsymbol{\Theta}$.

Next, we established similar asymptotic stationarity of RBMM for more than two
blocks. In this case, a well-known counterexample by Powell [50] shows that
block coordinates descent methods for $m\geq 3$ with block optimization may
not converge to stationary points. In deed, block coordinate descent is a
special case of block MM with the surrogate functions $g_{n}^{(i)}$ identical
to $f_{n}^{(i)}$. To overcome this issue, the proximal regularized version of
block coordinate descent was proposed ([31, 59, 2], which had the convergence
guarantees under certain assumptions. With this insight, we need to require
one more assumption on the surrogate function $g_{n}^{(i)}$. Namely, writing

$\displaystyle g_{n}^{(i)}(\theta)=f_{n}^{(i)}(\theta)+h_{n}^{(i)}(\theta),$
(22)

one can think of the nonnegative ‘surrogate gap’ function $h_{n}^{(i)}$ as a
regularization term. A primary example of such regularization in the Euclidean
case is the proximal point regularization, where we would take
$h_{n}^{(i)}(\theta)=\lambda\lVert\theta-\theta_{n-1}^{(i)}\rVert^{2}$. This
has the effect of penalizing a large change in the iterates
$\theta_{n}^{(i)}-\theta_{n-1}^{(i)}$. In the general Riemannian case, the
following assumption states that the surrogate gap function $h_{n}^{(i)}$ is
penalizing a large geodesic distance $d(\theta,\theta_{n-1}^{(i)})$.

######  (A1)iii.

(Distance-regularizing surrogates) There exists a strictly increasing function
$\phi:[0,\infty)\rightarrow\mathbb{R}$ such that $\phi(0)=0$ and

$\displaystyle
h_{n}^{(i)}(\theta):=g_{n}^{(i)}(\theta)-f_{n}^{(i)}(\theta)\geq
c\phi(d(\theta,\theta^{(i)}_{n-1}))$ (23)

for all $n\geq 1$ and $i=1,\dots,m$.

Note that if we use Riemannian proximal surrogates as in (A1)(i-rp), then
(A1)(iii) is automatically satisfied with $c=\lambda/2$ and $\phi(x)=x^{2}$.
This also holds with a possibly different parameter $c$ when we use the
Euclidean proximal surrogates within a compact constraint set on an embedded
manifold as in (A1)(i-ep). In fact, the geodesic distance and Euclidean
distance are equivalent over compact sets, see Lemma 4.10 for details. Also,
if the surrogates $g_{n}^{(i)}$ are $\rho$-strongly $g$-convex and if
$\theta_{n+1}^{(i)}$ is the exact minimizer of $g_{n}^{(i)}$, then by the
second-order growth property, one can verify

$\displaystyle
h_{n}^{(i)}(\theta)\geq\frac{\rho}{2}d^{2}(\theta,\theta^{(i)}_{n-1}).$ (24)

Hence (A1)(iii) is verified with $c=\rho/2$ and $\phi(x)=x^{2}$ in this case
as well.

Note a direct implication of (A1)(iii) is that if
$g_{n}^{(i)}(\theta_{n})-f_{n}^{(i)}(\theta_{n})=o(1)$ then
$d(\theta^{(i)}_{n},\theta^{(i)}_{n-1})=o(1)$. This fact will be crucial to
our proof of asymptotic stationarity of RBMM for $m\geq 3$ blocks, which is
our second main result stated below. Note that in the following result, we do
not require geodesic convexity of the constraint sets.

######  Theorem 3.3 (Asymptotic convergence to stationary points; many
blocks).

Let $f$ denote the objective function in (1) with $m\geq 1$. Let
$(\boldsymbol{\theta}_{n})_{n\geq 0}$ be an output of Algorithm 1. Under (A0),
(A1)(ii), and any of (i-gs), (i-rp) and (i-ep), every limit point of
$(\boldsymbol{\theta}_{n})_{n\geq 0}$ is a stationary point of $f$ over
$\boldsymbol{\Theta}$.

######  Remark 3.4 (Convexity of the constraint set).

The proof of Theorem 3.2 requires the constraint set
$\Theta^{(i)}\subseteq\mathcal{M}^{(i)}$ to be strongly convex w.r.t. the
geometry of the ambient manifold $\mathcal{M}^{(i)}$ (due to the Riemmanian
line search used in Prop. 6.13). However, Theorem 3.3 does not require the
constraint set $\Theta^{(i)}\subseteq\mathcal{M}^{(i)}$ to be strongly convex,
since we can use Assumption (A1)(iii) to avoid using Prop. 6.13. Therefore,
Theorem 3.3 applies when $\Theta^{(i)}$ is a manifold by itself (e.g., low-
rank manifold or Stiefel manifold) and $\mathcal{M}^{(i)}$ is a Euclidean
space in which $\Theta^{(i)}$ is embedded. See Section 4.2 for details.

Now we state our result concerning the rate of convergence of RBMM for the
case of Riemannian/Euclidean proximal (and not necessarily $g$-smooth)
surrogates.

######  Theorem 3.5 (Rate of convergence for Riemannian/Euclidean proximal
surrogates).

Let $f$ denote the objective function in (1) with $m\geq 1$. Let
$(\boldsymbol{\theta}_{n})_{n\geq 0}$ be an output of Algorithm 1 under (A0)
and (A1)(ii), and either (i-rp) or (i-ep). Assume the geodesic convexity of
the constraint sets. Then the following hold:

(i)

(Rate of convergence) There exists constants $M$ and $c>0$ independent of
$\boldsymbol{\theta}_{0}$ such that

$\min_{1\leq k\leq n}\left[-\sum_{i=1}^{m}\inf_{\eta\in
T_{\theta_{k}^{(i)}}^{*},\|\eta\|\leq
1}\left\langle\operatorname{grad}_{i}f(\boldsymbol{\theta}_{k}),\frac{\eta}{\min\\{r_{0},1\\}}\right\rangle\right]\leq\frac{M+c\sum_{n=1}^{\infty}\Delta_{n}}{\sqrt{n}/\log
n}$ (25)

(ii)

(Worst-case iteration complexity) The worst-case iteration complexity
$N_{\epsilon}$ for Algorithm 1 satisfies
$N_{\epsilon}=O(\varepsilon^{-2}\left(\log\varepsilon^{-2}\right)^{2})$

Theorem 3.5 establishes that RBMM with either Riemannian or Euclidean proximal
regularizer achieves an iteration complexity of
$\widetilde{O}(\varepsilon^{-2})$. The case of Euclidean proximal regularizer
in (A1)(i-ep) may be of wider practical interest than the Riemannian proximal
regularizer. Note in (A1)(i-rp), the Riemannian proximal regularizer is
geodesic distance squared, which may introduce additional computational
difficulties, especially when the close-form expression of end-point geodesic
distance is unknown. When the underlying manifold is embedded in Euclidean
space and the constraint set is compact (or the manifold itself is compact),
as in (A1)(i-ep), we can replace the Riemannian proximal regularizer by the
Euclidean proximal regularizer, which is often much easier to deal with
computationally.

Next, the theorem below states similar rates of convergence results for
$g$-smooth surrogates on general manifolds.

######  Theorem 3.6 (Rate of convergence for $g$-smooth surrogates).

Let $f$ denote the objective function in (1) with $m\geq 1$. Let
$(\boldsymbol{\theta}_{n})_{n\geq 0}$ be an output of Algorithm 1 under (A0)
and (A1)(ii), and (i-gs). Assume the geodesic convexity of the constraint
sets. Further assume that (A1)(iii) holds with $\phi(x)=cx^{2}$ for some
constant $c>0$. Then the following hold:

(i)

(Worst-case rate of convergence) There exists constants $M,c>0$ independent of
$\boldsymbol{\theta}_{0}$ such that

$\min_{1\leq k\leq n}\left[-\sum_{i=1}^{m}\inf_{\eta\in
T_{\theta_{k}^{(i)}}^{*},\|\eta\|\leq
1}\left\langle\operatorname{grad}_{i}f(\boldsymbol{\theta}_{k}),\frac{\eta}{\min\\{r_{0},1\\}}\right\rangle\right]\leq\frac{M+c\sum_{n=1}^{\infty}\Delta_{n}}{n^{1/4}/(\log
n)^{1/2}}$ (26)

(ii)

(Worst-case iteration complexity) The worst-case iteration complexity
$N_{\epsilon}$ for Algorithm 1 satisfies
$N_{\epsilon}=O\left(\varepsilon^{-4}\left(\log\varepsilon^{-2}\right)\right)$.

(iii)

(Optimal convergence rate) Further assume that the surrogate gaps
$h_{n}^{(i)}=g_{n}^{(i)}-f_{n}^{(i)}$ satisfy $h_{n}^{(i)}(\theta)\leq
Cd^{2}(\theta,\theta_{n}^{(i)})$ for some constant $C>0$. Then

$\min_{1\leq k\leq n}\left[-\sum_{i=1}^{m}\inf_{\eta\in
T_{\theta_{k}^{(i)}}^{*},\|\eta\|\leq
1}\left\langle\operatorname{grad}_{i}f(\boldsymbol{\theta}_{k}),\frac{\eta}{\min\\{r_{0},1\\}}\right\rangle\right]\leq\frac{M+c\sum_{n=1}^{\infty}\Delta_{n}}{n^{1/2}/(\log
n)^{1/2}}$ (27)

and the worst-case iteration complexity $N_{\epsilon}$ for Algorithm 1 is
$N_{\epsilon}=O\left(\varepsilon^{-2}\left(\log\varepsilon^{-2}\right)\right)$.

Lastly, we state a practical corollary of Theorem 3.6 for Riemannian
optimization problems involving Stiefel manifolds (see Ex. 4.1). A special
case of it was stated in Corollary 1.1 in the introduction. An important fact
about Stiefel manifolds is that, if a function $f$ is (Euclidean) $L$-smooth
in the ambient Euclidean space, then it is $g$-smooth with respect to the
geometry on the Stiefel manifold for some smoothness parameter $L^{\prime}$
(see Lem. 4.9). Therefore, when the underlying manifolds are either Euclidean
or Stiefel, we can apply Theorem 3.6 for Euclidean smooth objectives and
surrogates and obtain iteration complexity of
$\widetilde{O}(\varepsilon^{-2})$. This expands the applicability of our
result to various optimization problems involving Stiefel manifolds.

######  Corollary 3.7 (Complexity of RBMM on Stiefel manifolds).

Suppose each underlying manifold $\mathcal{M}^{(i)}$ is either a Stiefel
manifold or a Euclidean space. Assume the objective function $f$ in (1) is
Euclidean $L$-smooth for some $L>0$. Suppose the surrogates $g_{n}^{(i)}$ are
Euclidean $L^{\prime}$-smooth for some constants $L^{\prime}>0$ and for some
constant $c>0$,

$\displaystyle
h_{n}^{(i)}(\theta):=g_{n}^{(i)}(\theta)-f_{n}^{(i)}(\theta)\geq
c\lVert\theta-\theta_{n-1}^{(i)}\rVert^{2}$ (28)

for all $n\geq 1$ and $i=1,\dots,m$. Assume the constraint sets
$\Theta^{(1)},\dots,\Theta^{(m)}$ are geodesic convex. Allow inexact
computation in the sense of (A0)(ii) with $d(\cdot,\cdot)$ in (19) replaced by
the Euclidean distance. Then the iterates produced by Algorithm 1
asymptotically converge to the set of stationary points and the algorithm has
iteration complexity of $\widetilde{O}(\varepsilon^{-2})$.

It is important to note that the conditions one needs to check to apply
Corollary 3.7 are completely Euclidean, except the $g$-convexity of constraint
sets of Stiefel manifolds, which becomes vacuous when there are no additional
constraints on the Stiefel manifolds. The results presented in Corollary 3.7
are applicable to various MM methods on Stiefel manifolds, as discussed later
in Section 4.4, including the recent MM methods investigated in [16].
Furthermore, in Section 5.1 we study the geodesically constrained subspace
tracking problem as a stylized application.

##  4\. Applications

In this section, we discuss some examples of our general framework of RBMM and
its connection to other classical algorithms.

###  4.1. Examples of manifolds

In this section, we list several examples of manifolds that are typically used
in various machine learning problems.

######  Example 4.1 (Stiefel Manifold).

The Stiefel manifold $\mathcal{V}^{n\times k}$ is the set of all orthonormal
$k$-frames in $\mathbb{R}^{n}$. That is, it is the set of ordered orthonormal
$k$-tuples of vectors in $\mathbb{R}^{n}$, i.e.

$\mathcal{V}^{n\times k}=\left\\{A\in\mathbb{R}^{n\times
k}:A^{T}A=I_{k}\right\\},$ (29)

where $I_{k}$ denotes the $k\times k$ identity matrix. For
$X\in\mathbb{R}^{n\times k}$, denotes its SVD as $X=U\Sigma V^{T}$, then the
projection of $X$ onto $\mathcal{V}^{n\times k}$ is

$\textup{Proj}_{\mathcal{V}^{n\times k}}(X)=UV^{T}.$ (30)

$\blacktriangle$

######  Example 4.2 (Fixed-rank Matrices Manifold).

The set of matrices with fixed rank-$r$

$\mathcal{R}_{r}=\left\\{X\in\mathbb{R}^{m\times
n}:\operatorname{rank}(X)=r\right\\}$ (31)

is a smooth submanifold of $\mathbb{R}^{m\times n}$. For
$X\in\mathbb{R}^{m\times n}$, denote its SVD as $X=U\Sigma V^{T}$, where the
two matrices $U=\left[u_{1},u_{2},\ldots,u_{m}\right]$ and
$V=\left[v_{1},v_{2},\ldots,v_{n}\right]$ are orthogonal matrices. The
diagonal entries of $\Sigma$, which are the singular values of $X$, are
written in nonincreasing order,

$\sigma_{1}(X)\geq\sigma_{2}(X)\geq\cdots\geq\sigma_{\min\\{n,m\\}}(X)\geq 0.$
(32)

Then the projection of $X$ onto $\mathcal{R}_{r}$ is given by

$\textup{Proj}_{\mathcal{R}_{r}}(X)=\sum_{i=1}^{r}\sigma_{i}(X)u_{i}v_{i}^{T}:=U\Sigma_{r}V^{T}.$
(33)

$\blacktriangle$

Also, Hadamard manifolds are a class of manifolds that is widely studied in
the literature, since it includes many commonly encountered manifolds.
Hadamard manifolds are Riemannian manifolds with nonpositive sectional
curvature that are complete and simply connected, see [12] and [15]. Hadamard
manifolds have infinite injectivity radii at every point.

Below we provide some examples of Hadamard manifolds (more details can be
found in e.g. [9]).

######  Example 4.3 (Euclidean spaces).

The Euclidean space $\mathbb{R}^{n}$ with its usual metric is a Hadamard
manifold with constant sectional curvature equal to 0. $\blacktriangle$

######  Example 4.4 (Hyperbolic spaces).

We equip $\mathbb{R}^{n+1}$ with the $(-1,n)$-inner product

$\langle x,y\rangle_{(-1,n)}:=-x^{0}y^{0}+\sum_{i=1}^{n}x^{i}y^{i}$

for $x:=\left(x^{0},x^{1},\ldots,x^{n}\right)$ and
$y:=\left(y^{0},y^{1},\ldots,y^{n}\right)$. Define

$\mathbb{H}^{n}:=\left\\{x\in\mathbb{R}^{n+1}:\langle
x,x\rangle_{(-1,n)}=-1,x_{0}>0\right\\}.$

Then $\langle\cdot,\cdot\rangle$ induces a Riemannian metric $g$ on the
tangent spaces $T_{p}\mathbb{H}^{n}\subset T_{p}\mathbb{R}^{n+1}$ for
$p\in\mathbb{H}^{n}$. The sectional curvature of
$\left(\mathbb{H}^{n},g\right)$ is $-1$ at every point. $\blacktriangle$

######  Example 4.5 (Manifolds of positive definite matrices).

The space $\mathbb{S}_{++}^{n}$ of symmetric positive definite matrices
$n\times n$ with real entries is a Hadamard manifold if it is equipped with
the Riemannian metric

$\left\langle\Omega_{1},\Omega_{2}\right\rangle_{\Sigma}\triangleq\frac{1}{2}\operatorname{Tr}\left(\Omega_{1}\Sigma^{-1}\Omega_{2}\Sigma^{-1}\right)\quad\forall\Omega_{1},\Omega_{2}\in
T_{\Sigma}\mathbb{S}_{++}^{n}.$ (34)

$\blacktriangle$

###  4.2. Euclidean block MM

When specialized on the standard Euclidean manifold, our RBMM becomes the
standard Euclidean Block MM (e.g., see BSUM in [33]), where convergence rate
for convex and strongly convex objectives are known. Recently, Lyu and Li [39,
Thm. 2.1] obtained convergence rates for Euclidean Block MM algorithms for
nonconvex objectives with convex constraints. Our general results can recover
part of their results:

######  Corollary 4.6 (Complexity of Euclidean Block MM).

Theorems 3.2, 3.3, 3.5, and 3.6 hold for Algorithm 1 when the underlying
manifolds are Euclidean. In particular, the complexity result in Theorems 3.5
and 3.6 hold for the BSUM algorithm in [33].

We remark that [39, Thm. 2.1]] also covers the case when convex surrogates
with non-strongly-convex majorization gaps are used along with trust-regions
and diminishing radii. Corollary 4.6 does not cover this case.

Below we give some examples of the Euclidean block MM. One primary example is
Euclidean block proximal updates, namely, applying the surrogates in (20) when
the underlying manifold is the Euclidean space. The Euclidean block proximal
updates reads

$\theta_{n}^{(i)}\leftarrow\operatorname*{arg\,min}_{\theta\in\Theta^{(i)}}\left(g^{(i)}_{n}(\theta)=f^{(i)}_{n}(\theta)+\lambda_{n}\|\theta-\theta^{(i)}_{n-1}\|^{2}_{F}\right).$
(35)

Another example of Euclidean block MM is the following block prox-linear
update proposed in [59]: For minimizing a differentiable function $f$ defined
on the Euclidean space,

$\displaystyle\theta_{n}^{(i)}$
$\displaystyle\leftarrow\operatorname*{arg\,min}_{\theta\in\Theta^{(i)}}\left(g_{n}^{(i)}(\theta):=f_{n}^{(i)}(\theta_{n-1}^{(i)})+\langle\nabla
f_{n}^{(i)}(\theta_{n-1}^{(i)}),\,\theta-\theta_{n-1}^{(i)}\rangle+\frac{\lambda}{2}\lVert\theta-\theta_{n-1}^{(i)}\rVert^{2}\right).$
(36)

In [59], under a mild condition, it was shown that the above algorithm
converges asymptotically to a Nash equilibrium (a weaker notion than
stationary points) and also a local rate of convergence under the Kurdyka-
Łojasiewicz condition is established. Notice that when $f$ is block-wise
$L$-smooth and if $\lambda\geq L$, then $g_{n}^{(i)}$ is a majorizing
surrogate of $f_{n}^{(i)}$ at $\theta_{n-1}^{(i)}$. Thus (36) is a special
instance of our RBMM algorithm in this case and hence Theorems 3.2 and 3.6
apply. That is, our general results imply that the block prox-linear algorithm
(36) in the Euclidean space converges asymptotically to the stationary points
(not only Nash equilibrium) and also has iteration complexity of
$\widetilde{O}(\varepsilon^{-2})$.

In fact, the block prox-linear update (36) coincides with block projected
gradient descent with a fixed step size. Indeed, denoting $\nabla:=\nabla
f_{n}^{(i)}(\theta_{n-1}^{(i)})$, (36) is equivalent to

$\displaystyle\theta_{n}^{(i)}\leftarrow\operatorname*{arg\,min}_{\theta\in\Theta^{(i)}}\left(\langle\nabla,\,\theta\rangle+\frac{\lambda}{2}\lVert\theta\rVert^{2}-\lambda\langle\theta,\theta^{(i)}_{n-1}\rangle\right)$
$\displaystyle=\operatorname*{arg\,min}_{\theta\in\Theta^{(i)}}\left\lVert\theta-\left(\theta^{(i)}_{n-1}-\frac{1}{\lambda}\nabla\right)\right\rVert^{2}$
(37)
$\displaystyle=\textup{Proj}_{\Theta^{(i)}}\left(\theta^{(i)}_{n-1}-\frac{1}{\lambda}\nabla\right).$
(38)

Notice that since $\lambda\geq L$, the above becomes the standard block
projected gradient descent (Block PGD) update with step-size $\in(0,1/L]$. For
block PGD with convex objectives, convergence of function value with
complexity $\tilde{O}(\varepsilon^{-1})$ is established in [19]. Recently, in
[39], the authors showed the block PGD for smooth non-convex objectives
converges to the set of stationary points with complexity
$\tilde{O}(\varepsilon^{-2})$. Our general results apply to this classical
algorithm as well and recover the same complexity as [39].

When the constraint set $\Theta^{(i)}$ is a Riemannian manifold embedded in
the Euclidean space, in many cases it is non-convex as a subset of Euclidean
space. Therefore, when applying the prox-linear updates as in (36) and (37),
the complexity result in Theorem 3.5 and 3.6 do not hold, while the asymptotic
convergence result in Theorem 3.3 still hold. Moreover, The projection in (37)
can be solved exactly on some well-known manifolds, including low-rank
manifolds and Stiefel manifolds, which can be found in Section 4.1. More
examples and details can be found in [6].

Lastly, we provide two stylized applications of the Euclidean block MM in
Section 5, namely the robust PCA (Section 5.4) and the Riemannian CP-
dictionary learning (Section 5.3).

###  4.3. Block Riemannian proximal updates on Hadamard manifolds

We consider the Riemannian proximal surrogates in (A1)(i-rp) on Hadamard
manifolds. As shown later in Section 6.2, this Riemannian proximal surrogate
is geodesically strongly convex on Hadamard manifolds. Hence, the subproblem
in the minimization step in (LABEL:eq:RBMM_MmMm), i.e. minimizing
$g_{n}^{(i)}$ could be solved efficiently using classical Riemannian
optimization methods ([58], [64], [41]).

Our general results in Theorems 3.3 and 3.5 (with Riemannian proximal
surrogates) apply to Hadamard manifolds and we obtain the following corollary.

######  Corollary 4.7 (Complexity of block Riemannian proximal updates on
Hadamard manifolds).

Theorems 3.3 and 3.5 hold for block Riemannian proximal updates (20) on
Hadamard manifolds. That is, the RBMM with proximal surrogates converges
asymptotically to the set of stationary points and has iteration complexity of
$\widetilde{O}(\varepsilon^{-2})$.

In Section 5.2, we give a stylized example of block Riemannian proximal
updates solving the optimistic likelihood problem, where the manifold of PSD
matrices is involved.

###  4.4. Block Riemannian/Euclidean proximal and prox-linear updates on
Stiefel manifolds

In this section, we discuss both the Riemannian and Euclidean proximal/prox-
linear updates on the Stiefel manifold, as well as their variants and
applications.

It is known that the Stiefel manifold has non-negative sectional curvature
([62]), so it is not a Hadamard manifold. Nevertheless, compact manifolds have
a positive injectivity radius ([22, Thm. III.2.3]). In particular, the Stiefel
manifold has an injectivity radius of at least $0.89\pi$ ([52, Eq. (5.13)]),
which satisfies (A1)(ii). Hence the results in Theorems 3.3 and 3.5 apply and
we state the results in the following corollary.

######  Corollary 4.8 (Complexity of block Riemannian/Euclidean proximal
updates on Stiefel manifolds).

Theorems 3.3 and 3.5 hold for both block Riemannian proximal updates (20) and
block Euclidean proximal updates (21) on Stiefel manifolds. That is, the RBMM
with Euclidean/Riemannian proximal surrogates converges asymptotically to the
set of stationary points and has iteration complexity of
$\widetilde{O}(\varepsilon^{-2})$.

We remark that in Corollary 4.8, both Riemannian proximal surrogates and
Euclidean proximal surrogates are allowed. The Riemannian proximal surrogates
involve the geodesic distance in the updates, which brings additional
computational difficulty since the closed-form solution of the end-point
geodesic is unknown for the Stiefel manifold. A survey of numerical methods on
computing geodesics on the Stiefel manifold can be found in [27]. Instead, the
Euclidean proximal surrogates use the Euclidean distance function as a
regularizer, which provides computational savings.

In fact, our general framework of RBMM can utilize standard (Euclidean)
$L$-smooth surrogates for block optimization problems involving Stiefel
manifolds, as stated in Corollary 3.7. Below we first give a key lemma from
[21] for deriving Corollary 3.7 from Theorems 3.3 and 3.6. This lemma states
any $L$-smooth function in Euclidean space is a $g$-smooth function on the
Stiefel manifold, and therefore an $L$-smooth objective function and
surrogates satisfy the $g$-smoothness assumption in (A0) and (A1).

######  Lemma 4.9 (Euclidean smoothness implies geodesic smoothness on Stiefel
manifold; Lem. 2.4 and Appendix C.1 in [21]).

If $f$ is $L$-smooth in Euclidean space $\mathbb{R}^{n\times d}$, then there
exists a constant $L_{g}=L+L_{N}$ such that $f$ is $g$-smooth with parameter
$L_{g}$ on the Stiefel manifold $\mathcal{V}^{n\times d}$, where
$L_{N}=\max_{x\in\mathcal{V}^{n\times d}}\|\nabla f(x)\|$.

In order for the iteration complexity results in Corollary 3.7 to hold, one
also needs to verify assumption (A1)(iii). Namely, the majorization gap should
be lower bounded by an increasing function of geodesic distance. The following
geometric lemma states that for compact embedded submanifolds of the Euclidean
space, the Riemannian and the Euclidean distances are within constant
multiples of each other. Therefore, in this case, the quadratic majorization
gap in terms of the Euclidean distance in (28) implies the quadratic
majorization gap in terms of the Riemannian distance in (A1)(iii). This allows
us to use a wide range of computationally efficient surrogates on Stiefel
manifolds such as Euclidean proximal surrogates (21), Euclidean prox-linear
surrogates (36), and Euclidean regularized linear surrogates, which will be
discussed later in this section.

See Figure 3 for an illustration of Lemma 4.10 on Stiefel manifolds.

######  Lemma 4.10 (Equivalence of distance on Riemannian manifold; Lem. 4.1
in [45]).

Let $M\subset\mathbb{R}^{n}$ be a smooth submanifold, equipped with a
Riemannian metric $g$. The geodesic distance between $x,y\in M$ induced by $g$
is denoted as d(x,y). Consider the Euclidean norm $\|\cdot\|$ on
$\mathbb{R}^{n}$. Let $K\subseteq M$ be compact. Then there exists $c>0$ such
that for $x,y\in K$,

$cd(x,y)\leq\|x-y\|\leq d(x,y).$ (39)

Figure 3. Illustration of distances on $\mathcal{V}^{n\times p}$ with
$n=2,p=1$. The length of the blue segment is Euclidean distance $\|X-Y\|_{F}$.
The length of the arc between $X$ and $Y$ is $d(X,Y)$.

Below we give a proof of Corollary 3.7. Recall that Corollary 1.1 is a direct
consequence of Corollary 3.7.

###### Proof of Corollary 3.7.

To deduce Corollary 3.7 from Theorems 3.3 and 3.6, we verify the assumptions
one by one. First, by the hypothesis on Euclidean smoothness of the objective
$f$ and the surrogates $g_{n}^{(i)}$ and Lemma 4.9 (and also using the
definition of $d$ on the product manifold in (9)), we have that $f$ and
$g_{n}^{(i)}$s are also $g$-smooth. This verifies (A0)(i) and (A1)(i-gs). By
Lemma 4.10 and the hypothesis in Corollary 3.7, we verify (A0)(ii) and
(A1)(iii). Lastly, (A1)(ii) is satisfied since both Stiefel manifolds and
Euclidean space have a uniformly positive injectivity radius. Then the
assertion in Corollary 3.7 follows from Theorems 3.3 and 3.6. ∎

###  4.5. Iteration complexity of MM on Stiefel manifolds in [16]

In order to solve minimization over Stiefel manifolds

$\displaystyle\min_{U\in\mathcal{V}^{p\times k}}f(U),$ (40)

Breloy et al. [16] recently proposed an MM method incorporating linear
surrogates and reducing the surrogate minimization problem to a projection
onto the Stiefel manifold. More precisely, given an iterate
$U_{n-1}\in\mathcal{V}^{p\times k}$ at iteration $n-1$, the surrogate $g_{n}$
at iteration $n$ takes the following form

$g_{n}(U)=-\operatorname{Tr}\left(\mathbb{R}^{H}(U_{n-1})U\right)-\operatorname{Tr}\left(U^{H}\mathbb{R}(U_{n-1})\right)+\textup{const},$
(41)

where H is the conjugate transpose operator and
$\mathbb{R}(\cdot):\mathbb{C}^{p\times k}\to\mathbb{C}^{p\times k}$ is a
matrix function chosen so that $g_{n}(U_{n-1})=f(U_{n-1})$ and $g_{n}\geq f$.
Since $g_{n}$ is linear and since $U\in\mathcal{V}^{p\times k}$, we have

$\displaystyle
U_{n+1}\overset{\textup{def}}{=}\operatorname*{arg\,min}_{U\in\mathcal{V}^{p\times
k}}g_{n}(U)=\operatorname*{arg\,min}_{U\in\mathcal{V}^{p\times
k}}\,\,\lVert\mathbf{R}^{H}(U_{n-1})-U\rVert_{F}.$ (42)

When $\mathbf{R}^{H}(U_{n-1})$ is of full-rank, then the rightmost projection
problem has a unique solution given by the projection operator in (30). See
various applications of this MM method in [16].

In [16], the authors showed this MM method asymptotically converges to the set
of stationary points by adapting a convergence result of Euclidean BMM
algorithm in [53]. However, there is no known bound on the iteration
complexity due to the non-convexity of the constraint set and the objective
function. An iteration complexity of $\tilde{O}(\varepsilon^{-2})$ of
Euclidean BMM for non-convex smooth objectives with convex constraints has
been obtained recently in [39]. However, this result is not applicable here
since the Stiefel manifolds cannot be viewed as convex constraint sets within
Euclidean spaces.

By applying Corollary 3.7, we can obtain an iteration complexity bound of
$\tilde{O}(\varepsilon^{-2})$ for the above MM method on the Stiefel manifold
with a slight modification via Euclidean proximal regularization. Namely, for
a fixed proximal regularization parameter $\lambda\geq 0$, consider the
following surrogate

$\displaystyle\tilde{g}_{n}(U)$ $\displaystyle=g_{n}(U)+\lambda\lVert
U-U_{n-1}\rVert_{F}^{2}$ (43)
$\displaystyle=-\operatorname{Tr}\left((\mathbb{R}(U_{n-1})+\lambda
U_{n-1})^{H}U\right)-\operatorname{Tr}\left(U^{H}(\mathbb{R}(U_{n-1})+\lambda
U_{n-1}))\right)+\textup{const}.$ (44)

Note that minimizing the above surrogate $\tilde{g}_{n}$ is as easy as
minimizing $g_{n}$:

$\displaystyle
U_{n+1}\overset{\textup{def}}{=}\operatorname*{arg\,min}_{U\in\mathcal{V}^{p\times
k}}\tilde{g}_{n}(U)$
$\displaystyle=\operatorname*{arg\,min}_{U\in\mathcal{V}^{p\times
k}}\,\,\lVert\mathbf{R}^{H}(U_{n-1})+\lambda U_{n-1}-U\rVert_{F}.$ (45)

Note that the MM update in (45) satisfies the hypothesis of Corollary 3.7.
Indeed, the surrogates $g_{n}$ in (41) and $\tilde{g}_{n}$ in (43) are linear
and hence are $L$-smooth in Euclidean space. The addition of Euclidean
proximal regularization ensures that we have at least a quadratic majorization
gap as in (28). Therefore, by Corollary 3.7, the iteration complexity of the
MM update in (45) is of $\tilde{O}(\varepsilon^{-2})$. These results are
formally stated in Corollary 4.11.

######  Corollary 4.11 (Complexity of MM with regularized linear surrogates on
the Stiefel manifold).

Consider the problem of minimizing a differentiable objective $f$ on the
Stiefel manifold $\mathcal{V}^{p\times k}$ as in (40).

Then the iterates generated by (45) with arbitrary initialization converges
asymptotically to the set of stationary points of (40). Moreover, the
iteration complexity is $\widetilde{O}(\varepsilon^{-2})$.

To the best of our knowledge, this is the first complexity result for MM on
Stiefel manifolds in the literature. In particular, the iteration complexity
in Corollary 4.11 applies for various problem instances discussed in [16, Sec.
V] including power iteration for computing top eigenvector, generic non-
homogenous quadratic form, and nonconvex subspace recovery. In Section 5.1, we
use a similar idea to obtain the same iteration complexity bound for block MM
algorithm for the geodesically constrained subspace tracking problem [18].

##  5\. Stylized Applications

In this section, we will discuss the following four stylized applications of
our theory:

  1. (1)

Geodesically constrained subspace tracking [18];

  2. (2)

Optimistic likelihood under Fisher-Rao distance [49];

  3. (3)

Riemannian CP-dictionary-learning [42, 26];

  4. (4)

Robust PCA [23, 55].

The first problem above is an application of block Euclidean proximal updates
on Stiefel manifolds, which verifies our assumptions for Theorem 3.6 and we
are able to derive a new iteration complexity result (Cor. 5.1), as
aforementioned in Corollary 3.7 and Section 4.4. The second problem above is
an application of Riemannian proximal updates on Hadamard manifolds (see
Section 4.3), which verifies our assumptions for Theorem 3.5 and we are able
to derive a new iteration complexity result (Cor. 5.2). For the last two
problems above, we only prove asymptotic convergence to stationary points by
using Theorem 3.3. These two problems are formulated as minimizing a cost
function involving the Euclidean distance function (e.g., the matrix Frobenius
norm) over low-rank manifolds. The Euclidean distance function is not
$g$-smooth over low-rank manifolds (see Appendix C.1 for details). Therefore,
in order to satisfy the assumption of $g$-smoothness of the objective in
(A0)(i), we choose Euclidean geometry as the underlying manifold structure,
and take the embedded submanifolds as constraints. In this case, the
constraints are not $g$-convex with respect to the underlying Euclidean
geometry, so we are not able to apply our iteration complexity results
(Theorems 3.5 and 3.6). However, we can still deduce asymptotic convergence to
stationary points by using Theorem 3.3 (see Section 4.2), which fortunately
does not require $g$-convexity of the constraint sets (see Remark 3.4).

###  5.1. Geodesically constrained subspace tracking [18]

Let $X_{i}\in\mathbb{R}^{d\times l}$ for $i=1,\cdots,T$ be data generated from
a low-rank model with noise,

$X_{i}=U_{i}G_{i}+N_{i},$ (46)

where $U_{i}\in\mathbb{R}^{d\times k}$ has orthonormal columns representing a
point on the Grassmannian $\mathcal{G}(k,d)$, the space of all rank-$k$
subspaces in $\mathbb{R}^{d};G_{i}\in\mathbb{R}^{k\times\ell}$ holds weight or
loading vectors; and $N_{i}\in\mathbb{R}^{d\times\ell}$ is an independent
additive noise matrix.

For the geodesic subspace tracking problem, we observe $X_{i}$ and our
objective is to estimate $U_{i}$ for $i=1,\cdots,T$. Let $\mathcal{V}^{d\times
k}$ denote the Stiefel manifold. We model each $U_{i}$ as an orthonormal basis
whose span is sampled from a single continuous Grassmanian geodesic
$U(t):[0,1]\to\mathcal{V}^{d\times k}$, parameterized as follows: For
$H,Y\in\mathcal{V}^{d\times k}$ consisting of orthogonal columns, i.e.,
$H^{T}Y=O$,

$U_{i}=U(t_{i})=H\cos(\Theta t_{i})+Y\sin(\Theta t_{i}).$ (47)

Here $\Theta\in\mathbb{R}^{k\times k}$ is a diagonal matrix where its $j$th
diagonal entry, $\theta_{j}$, is the $j$-th principal angle between the two
endpoints of the geodesic. We assume either the time-points $t_{i}$ are given,
or the observed matrices $X_{i}$ are equidistant along a geodesic curve.

The objective function is formulated as follows,

$f(U)=f(H,Y,\Theta)=\min_{\\{G_{i}\\}_{i=1}^{T}}\|X_{i}-U(t_{i})G_{i}\|^{2}_{F}=-\sum_{i=1}^{T}\|X_{i}^{T}U(t_{i})\|^{2}_{F}+c,$
(48)

where for the last equality we have substituted the optimal
$G_{i}=U\left(t_{i}\right)^{T}X_{i}$ and $c$ is a constant (see [30]).

To approximately minimize (48), Blocker et al. proposed a two-block MM
approach in [18], where one alternatively optimizes two block parameters
$Q:=(H,Y)\in\mathcal{V}^{d\times 2k}$ and $\Theta$. Here we propose a proximal
regularized version of the BMM method in [18] and establish its asymptotic
convergence property in Corollary 5.1. For proximal regularization parameters
$\lambda_{Q},\lambda_{\Theta}\geq 0$, the proposed algorithm reads as

$\displaystyle Q_{n+1}\leftarrow WV^{T},\,\,\textup{where $W\Sigma V^{T}$ is
the SVD of $\lambda_{Q}Q_{n}-\nabla_{Q}f(Q_{n},\Theta_{n})$}$ (49)
$\displaystyle(\theta_{j})_{n+1}\leftarrow(\theta_{j})_{n}-\frac{1}{w_{j}((\theta_{j})_{n})+\lambda_{\Theta}}\nabla
f^{(2)}_{n+1,j}\left((\theta_{j})_{n}\right)\,\,\textup{for $j=1,\dots,k$},$
(50)

where $(\theta_{j})_{n}$ denotes the $j$th diagonal entry of the $k\times k$
diagonal matrix $\Theta_{n}$ at iteration $n$,
$w_{j}=\sum_{i=1}^{T}w_{f_{i,j}}$ is the "weighting function" defined in [18],
and

$f^{(2)}_{n+1,j}\left(\theta_{j}\right):=-\sum_{i=1}^{T}(r_{i,j})_{n+1}\cos\left(2\theta_{j}t_{i}-(\phi_{i,j})_{n+1}\right)+(b_{i,j})_{n+1}.$
(51)

The definition of the parameters $\phi_{i,j}$, $r_{i,j}$ and $b_{i,j}$ in (51)
can be found in [18], which we also provide in Appendix C.2 for completeness.
Note that when $\lambda_{Q}=0$, (49) becomes the updates for $Q$ in [18].
Similarly, $\lambda_{\Theta}=0$ gives the updates for $\Theta$ in [18]. Thus
our algorithm (49)-(50) generalizes the BMM algorithm proposed in [18].

We first argue why we can cast the above algorithm as RBMM. The discussion we
provide here is a minor modification of the derivation in [18]. First, we
discuss the $Q$-update (49). Let $Z_{i}=[\cos(\Theta t_{i})\;;\;\sin(\Theta
t_{i})]$, which is the vertical concatenation of $\cos(\Theta t_{i})$ and
$\sin(\Theta t_{i})$. The objective function can be rewritten as

$f(Q,\Theta)=-\sum_{i=1}^{T}\|X_{i}^{T}QZ_{i}\|_{F}^{2}$ (52)

and the gradient with respect to the first block $Q$ is given by
$\nabla_{Q}f=-2\sum_{i=1}^{T}X_{i}X_{i}^{T}QZ_{i}Z_{i}^{T}$. The marginal
objective function for updating $Q$ can be rewritten as

$\displaystyle f_{n+1}^{(1)}(Q)$
$\displaystyle:=-\left\langle\sum_{i=1}^{T}X_{i}X_{i}^{T}Q(Z_{i})_{n}(Z_{i})_{n}^{T},Q\right\rangle,$
(53)

which is a concave-up quadratic function in $Q$. Also note that
$f_{n+1}^{(1)}$ is $L$-smooth for some constant $L>0$ over the compact
parameter space.

We consider the following proximal majorizer of $f_{n+1}^{(1)}$: For
$\lambda\geq 0$,

$\displaystyle g_{n+1}^{(1)}(Q)$
$\displaystyle:=-\left\langle\sum_{i=1}^{T}X_{i}X_{i}^{T}Q_{n}(Z_{i})_{n}(Z_{i})_{n}^{T},Q\right\rangle+\frac{\lambda}{4}\lVert
Q-Q_{n}\rVert_{F}^{2}$ (54)
$\displaystyle=\frac{1}{2}\langle\nabla_{Q}f(Q_{n},\Theta_{n}),Q\rangle+\frac{\lambda}{4}\lVert
Q-Q_{n}\rVert_{F}^{2}.$ (55)

We claim that

$\displaystyle Q_{n+1}$
$\displaystyle=\operatorname*{arg\,min}_{Q\in\mathcal{V}^{d\times
2k}}g_{n+1}^{(1)}(Q)=\begin{cases}\textup{Proj}_{\mathcal{V}^{d\times
2k}}\left(Q_{n}-\frac{1}{\lambda}\nabla_{Q}f(Q_{n},\Theta_{n})\right)&\textup{if
$\lambda>0$}\\\ \textup{Proj}_{\mathcal{V}^{d\times
2k}}\left(-\nabla_{Q}f(Q_{n},\Theta_{n})\right)&\textup{if
$\lambda=0$}.\end{cases}$ (56)

Indeed, the above MM update for $\lambda>0$ follows from (37). For
$\lambda=0$, we use the fact that $Q^{T}Q=I$ for all $Q$ in the Stiefel
manifold so that

$\displaystyle\operatorname*{arg\,min}_{Q\in\mathcal{V}^{d\times
2k}}\langle\nabla_{Q}f(Q_{n},\Theta_{n}),Q\rangle=\operatorname*{arg\,min}_{Q\in\mathcal{V}^{d\times
2k}}\lVert
Q+\nabla_{Q}f(Q_{n},\Theta_{n})\rVert^{2}_{F}=\textup{Proj}_{\mathcal{V}^{d\times
2k}}\left(-\nabla_{Q}f(Q_{n},\Theta_{n})\right).$ (57)

Notice that projecting onto the Stiefel manifold can be easily done by SVD
(see, e.g., [6]). Hence (56) coincides with (49).

Next, we discuss the $\Theta$-update. For conciseness, we only put the
expression of the loss function and surrogates here. For full details, we
refer the readers to [18]. The marginal loss function for $\Theta$ is
separable for each diagonal element $\theta_{j}$ of $\Theta$. Consider the
following proximal majorizer of $f_{n+1,j}^{(2)}$ in (51): For
$\lambda_{\Theta}\geq 0$,

$g_{n+1,j}^{(2)}(\theta_{j})=f^{(2)}_{n+1,j}\left((\theta_{j})_{n}\right)+\nabla
f^{(2)}_{n+1,j}\left((\theta_{j})_{n}\right)\left(\theta_{j}-(\theta_{j})_{n}\right)+\frac{w_{j}((\theta_{j})_{n})+\lambda_{\Theta}}{2}\left(\theta_{j}-(\theta_{j})_{n}\right)^{2}$
(58)

where $(\theta_{j})_{n}$ is the value of $\theta_{j}$ at iteration $n$. Then
by using (37), we see that (50) coincides with minimizing the majorizing
surrogate $g_{n+1,j}^{(2)}$ above.

Now we discuss the convergence of this block MM algorithm. In fact, we can
apply Corollary 3.7 and directly get the convergence and complexity results.
The asymptotic convergence result in Corollary 3.7 for the updates (49)-(50)
with $\lambda_{Q},\lambda_{\Theta}\geq 0$. Note when
$\lambda_{Q}=\lambda_{\Theta}=0$, the updates (49)-(50) become the vanilla
block MM method in [18]. Moreover, when the proximal parameter
$\lambda_{Q},\lambda_{\Theta}$ are strictly positive, the complexity result in
Corollary 3.7 holds. We state the convergence and complexity results of this
proximal regularized MM method in the following corollary:

######  Corollary 5.1 (Convergence and complexity of regularized BMM for
geodeiscally constrained subspace tracking).

Given a sequence of data $X_{i}$ for $i=1,\cdots,T$. Let
$\boldsymbol{\theta}_{k}=(Q_{k},\Theta_{k})$ be generated by (49)-(50) with
arbitrary initialization
$\boldsymbol{\theta}_{0}\in\boldsymbol{\Theta}=\mathcal{V}^{d\times
2k}\times\mathbb{R}^{k\times k}$. Suppose the proximal parameters
$\lambda_{Q}$ and $\lambda_{\Theta}$ are non-negative. Then the limit points
of $(\boldsymbol{\theta}_{n})_{n\geq 0}$ are stationary points of problem
(52). Moreover, if $\lambda_{Q},\lambda_{\Theta}>0$, then the iteration
complexity is $\widetilde{O}(\varepsilon^{-2})$.

###### Proof.

In fact, in order to apply Corollary 3.7, we only need to verify the Euclidean
smoothness of the marginal loss function. Note

$\displaystyle\nabla_{Q}f$
$\displaystyle=-2\sum_{i=1}^{T}X_{i}X_{i}^{T}QZ_{i}Z_{i}^{T},\qquad\nabla
f^{(2)}_{n+1,j}\left(\theta_{j}\right)=2\sum_{i=1}^{T}(r_{i,j})_{n+1}t_{i}\sin\left(2\theta_{j}t_{i}-(\phi_{i,j})_{n+1}\right),$
(59)

so $f$ is block-wise Euclidean $L$-smooth. Hence, the convergence and
complexity results in Corollary 3.7 follow. ∎

Now we compare our results with other existing works. We remark that the BMM
algorithm in [18] is based on the MM algorithm on Stiefel manifold by Breloy
et al. [16]. As aforementioned in Section 4.4, the asymptotic convergence to
the set of all stationary points in [16] was established by adopting a
convergence result of Euclidean BMM algorithm in [53]. The complexity results
are still unknown due to the non-convexity of the constraint set. Here, we
give the first complexity result in the literature. Moreover, our method is
computationally efficient with close-form updates as shown in (49)-(50).

(a)

(b)

(c)

(d)

Figure 4. Convergence of RBMM in geodesic error under different settings.
Average geodesic error is computed over 50 independent trials. The dimension
is $d=30$ and the additive Gaussian noise has a standard deviation
$\sigma=0.1$. The value of other parameters is shown in the title for each
panel. Adding Euclidean proximal regularizer does not accelerate the
convergence, mainly due to the orthogonal property on Stiefel manifolds.
Namely, any point $Q$ on the Stiefel manifold satisfies $Q^{T}Q=\mathbf{I}$.
Hence the proximal regularizer degenerates into a linear term providing with
limited acceleration.

Next, we discuss the results of numerical experiments. In the case of
synthetic data, we have access to the true geodesic, so we could compare the
estimated geodesic $\hat{U}(t)$ with the true geodesic $U(t)$. The following
error metric is used, which is the square root of the average squared subspace
error between corresponding points along the geodesic,

$\text{ Geodesic Error
}=\sqrt{\int_{0}^{1}\frac{1}{2k}\left\|\hat{U}(t)\hat{U}(t)^{T}-U(t)U(t)^{T}\right\|_{\mathrm{F}}^{2}\mathrm{\leavevmode\nobreak\
d}t}.$ (60)

In practice, we approximate the integral by Riemann sum using the sample
points. The geodesic error (60) takes minimum value of 0 when
$\operatorname{span}(\hat{U})=\operatorname{span}(U)$ and maximum value of 1
when $\operatorname{span}(\hat{U})\perp\operatorname{span}(U)$.

In the numerical experiments, we set the dimension $d=30$ and set the elements
of the noise matrix $N_{i}$ to be independent Gaussian noise with standard
deviation $\sigma=0.1$. The other parameters including the rank $k$, number of
samples $l$, number of time points $T$, and proximal parameters $\lambda_{Q}$
and $\lambda_{\Theta}$ are set differently for different experiments. We run
each experiment 50 times independently to compute the average geodesic error.
As shown in Figure 4, RBMM with different $\lambda$ converges under different
settings and is as good in terms of convergence speed as the block MM method
in [18]. Furthermore, Figure 4 shows that changing the proximal parameter
$\lambda_{Q}$ and $\lambda_{\Theta}$ only slightly affects the rate of
convergence, this is due to the following two reasons: first, there are only
two blocks for this problem, which gives a relatively simple setting for RBMM;
second, note $Q\in\mathcal{V}^{d\times 2k}$ which gives $Q^{T}Q=I_{2k}$,
therefore the quadratic term added would reduce to a linear term, which does
not significantly accelerate convergence. For more numerical results of this
geodesic subspace tracking problem, we refer the reader to [18].

###  5.2. Optimistic likelihood under Fisher-Rao distance [49]

Given a set of i.i.d data points $x^{M}\triangleq
x_{1},\ldots,x_{M}\in\mathbb{R}^{N}$ that are generated from one of several
Gaussian distributions $\mathbb{P}_{c}$, where $c\in C$ and $|C|<\infty$. We
want to determine the distribution $\mathbb{P}_{c^{*}}$, $c^{*}\in C$, under
which the following log-likelihood function $\ell(x^{M},\mathbb{P}_{c})$ is
maximized, i.e. we want to solve

$c^{*}=\operatorname*{arg\,max}_{c\in
C}\left\\{\ell\left(x^{M},\mathbb{P}_{c}\right)\triangleq-\frac{1}{M}\sum_{m=1}^{M}\left(x_{m}-\mu_{c}\right)^{T}\Sigma_{c}^{-1}\left(x_{m}-\mu_{c}\right)-\log\operatorname{det}\Sigma_{c}\right\\},$
(61)

where $\mu_{c}$ and $\Sigma_{c}$ denote the means and covariance matrices of
$\mathbb{P}_{c}$. Since methods that sample points to get good estimators of
$\mu_{c}$ and $\Sigma_{c}$ are usually costly, we consider the following
optimistic likelihood problem instead. Namely, instead of aiming to find the
Gaussian distribution from the set, we look for a Gaussian distribution that
maximizes the likelihood close to the empirical distribution under some
distance measure. More precisely, we consider the following optimistic
likelihood problem

$\max_{\mathbb{P}\in\mathcal{P}_{c}}\ell\left(x^{M},\mathbb{P}\right)\text{
with
}\mathcal{P}_{c}=\left\\{\mathbb{P}\in\mathcal{P}:\varphi\left(\hat{\mathbb{P}}_{c},\mathbb{P}\right)\leq\rho_{c}\right\\},$
(62)

where $\mathcal{P}$ is the set of all non-degenerate Gaussian distributions on
$\mathbb{R}^{N}$, $\hat{\mathbb{P}}_{c}$ is the empirical distribution
estimated from training data, $\varphi$ is the Fisher-Rao distance, and
$\rho_{c}\in\mathbb{R}_{+}$are the radii of the ambiguity sets
$\mathcal{P}_{c}$. For conciseness, we put the details of the Fisher-Rao
distance in Appendix C.3. We refer the readers to [49] for full details.

Now denoting the empirical mean and covariance from the data by $\hat{\mu}$
and $\hat{\Sigma}$, we can explicitly write down the optimization problem as

$\min_{\mu,\Sigma}f(\mu,\Sigma)\triangleq\left\langle
M^{-1}\sum_{m=1}^{M}\left(x_{m}-\mu\right)\left(x_{m}-\mu\right)^{T},\Sigma^{-1}\right\rangle+\log\operatorname{det}\Sigma,$
(63)

where
$\mu\in\Theta^{(1)}=\left\\{\mu\in\mathbb{R}^{N}:(\mu-\hat{\mu})^{T}(\mu-\hat{\mu})\leq\rho_{1}^{2}\right\\}$
and
$\Sigma\in\Theta^{(2)}=\left\\{\Sigma\in\mathbb{S}_{++}^{N}:d(\Sigma,\hat{\Sigma})\leq\rho_{2}\right\\}$.
Here $d(\Sigma,\hat{\Sigma})$ is the Fisher-Rao distance of two Gaussian
distributions with identical mean (see Appendix C.3).

Denote $\mu_{n}$ and $\Sigma_{n}$ as the $\theta_{n}^{(1)}$ and
$\theta_{n}^{(2)}$ generated by Algorithm 1 applied on $f(\mu,\Sigma)$. The
marginal objective functions are denoted by

$\displaystyle f^{(1)}_{n}:=f(\mu,\Sigma_{n-1})=\left\langle
M^{-1}\sum_{m=1}^{M}\left(x_{m}-\mu\right)\left(x_{m}-\mu\right)^{T},\Sigma^{-1}_{n-1}\right\rangle+\log\operatorname{det}\Sigma_{n-1}$
(64) $\displaystyle f^{(2)}_{n}:=f(\mu_{n},\Sigma)=\left\langle
S_{n},\Sigma^{-1}\right\rangle+\log\operatorname{det}\Sigma\qquad\text{where}\qquad
S_{n}=M^{-1}\sum_{m=1}^{M}\left(x_{m}-\mu_{n}\right)\left(x_{m}-\mu_{n}\right)^{T}.$
(65)

As aforementioned in Example 4.5, this manifold of positive definite matrices
is a Hadamard manifold. Thus we could construct Riemannian proximal surrogate
functions as in (20), more discussions can be found in Section 6.2, i.e.

$\displaystyle\mu_{n}=\operatorname*{arg\,min}_{\mu\in\mathbb{R}^{N}}g_{n}^{(1)}(\mu):=\left\langle
M^{-1}\sum_{m=1}^{M}\left(x_{m}-\mu\right)\left(x_{m}-\mu\right)^{T},\Sigma^{-1}_{n-1}\right\rangle+\log\operatorname{det}\Sigma_{n-1}+\frac{\lambda}{2}\|\mu-\mu_{n-1}\|^{2}$
(66)
$\displaystyle\Sigma_{n}=\operatorname*{arg\,min}_{\Sigma\in\mathbb{S}^{N}_{++}}g_{n}^{(2)}(\Sigma):=\left\langle
S_{n},\Sigma^{-1}\right\rangle+\log\operatorname{det}\Sigma+\frac{\lambda}{4}\left\|\log\left(\Sigma_{n-1}^{-\frac{1}{2}}\Sigma\Sigma_{n-1}^{-\frac{1}{2}}\right)\right\|_{F}^{2}.$
(67)

Note this block MM is a special instance of proximal methods on the Hadamard
manifold discussed in Section 4.3. Hence Theorems 3.3 and 3.5 apply. We state
it as a corollary as follows,

######  Corollary 5.2 (Complexity of Riemannian proximal updates for
optimistic likelihood problem).

Given a set of data points $(x_{i})_{i}$ for $i=1,\cdots,M$. Let
$\boldsymbol{\theta}_{n}=(\mu_{n},\Sigma_{n})$ be the iterates generated by
(66) with arbitrary initialization
$\boldsymbol{\theta}_{0}\in\boldsymbol{\Theta}\subseteq\mathbb{R}^{N}\times\mathbb{S}^{N}_{++}$.
Suppose the proximal regularization parameter $\lambda$ is strictly positive.
Then the limit points of $(\boldsymbol{\theta}_{n})_{n}$ are stationary points
of problem (63). Furthermore, it has iteration complexity of
$\widetilde{O}(\varepsilon^{-2})$.

###### Proof.

Note we only need to show (A0), (A1), and also the geodesic convexity of the
constraint set are satisfied. The proof of (A0)(i) and geodesic convexity of
the constraint sets are established based on some propositions adapted from
[49] which can be found in Appendix C.3. (A1)(i-rp) and (A0)(ii) holds by our
choice of Riemannian proximal surrogates. (A1)(ii) holds since the underlying
manifolds are Hadamard manifolds. Hence the corollary holds as a result of
Theorems 3.3 and 3.5. ∎

(a)

(b)

Figure 5. Comparison of block minimization and RBMM applied to optimistic
likelihood problem under Fisher-Rao distance. RBMM is implemented with
$\lambda=0.01,0.1,1$ respectively. Panel (A) shows the averaged error for the
first 300 iterations and Panel (B) shows the averaged error using a log-log
plot for 1000 iterations.

Now we show some numerical results. For numerical experiments, we use the same
setup as in [49] and study the empirical convergence behavior. We compare the
performance of block minimization using $f_{k}^{(i)}(\theta)$ for $i$-th block
and RBMM using $g_{k}^{(i)}(\theta)=f_{k}^{(i)}(\theta)+\lambda
d^{2}(\theta,\theta_{k-1}^{(i)})$ for $i$-th block with different values of
$\lambda$. We set the dimension of the data to be $N=10$ and denote
$f_{k}=f(\mu_{k},\Sigma_{k})$, the relative improvement is thus denoted as
$|f_{k+1}-f_{k}|/f_{k+1}$, which is computed via 10 independent experiments.
Numerical results are shown in Figure 5. While RBMM and block minimization
both perform well, RBMM is slightly faster, where the fastest convergence of
RBMM is achieved with $\lambda=0.01$.

###  5.3. Riemannian CP-dictionary-learning

(a) Euclidean: Synthetic data 1

(b) Euclidean: Synthetic data 2

(c) Stiefel: Synthetic data 1

(d) Low-rank: Synthetic data 1

Figure 6. (A) and (B) are some typical cases of synthetic data in Euclidean
case; (C) is the typical result when the first block is Stiefel manifold; (D)
is a synthetic example where the first block is a point on low-rank manifold.
The average relative reconstruction error with standard deviations are shown
by the solid lines and shaded regions of respective colors.

In the CANDECOMP/PARAFAC (CP) decomposition problem [36], given a data tensor
$X\in\mathbb{R}^{I_{1}\times\dots\times I_{m}}$ and an integer $R>0$, we would
like to find the loading matrices $U^{(i)}\in\mathbb{R}^{I_{i}\times R}$ for
$i=1,\cdots,m$ such that

$X\approx\sum_{k=1}^{R}\bigotimes_{i=1}^{m}U^{(i)}[:,k],$ (68)

where $U^{(i)}[:,k]$ denotes the $k^{\text{th }}$ column of the $I_{i}\times
R$ loading matrix matrix $U^{(i)}$ and $\otimes$ denotes the outer product. We
could formulate the above tensor decomposition problem as the following
optimization problem:

$\underset{U^{(1)}\in\Theta^{(1)},\ldots,U^{(m)}\in\Theta^{(m)}}{\operatorname{argmin}}\left(f\left(U^{(1)},\ldots,U^{(m)}\right):=\left\|X-\sum_{k=1}^{R}\bigotimes_{i=1}^{m}U^{(i)}[:,k]\right\|_{F}^{2}\right),$
(69)

where $\Theta^{(i)}\subseteq\mathbb{R}^{I_{i}\times R}$ is an embedded
manifold, which gives Riemannian constraints on the factor matrices. This
setup of Riemmanian CP-dictionary learning is related to the recent work [26],
where the authors used a CP-decomposition with Riemannian structure on the
space of factor matrices as a pre-conditioning algorithm for tensor
completion.

It is easy to see (69) is equivalent to

$\underset{U^{(1)}\in\Theta^{(1)},\ldots,U^{(m)}\in\Theta^{(m)}}{\operatorname{argmin}}\left\|X-\operatorname{Out}\left(U^{(1)},\ldots,U^{(m-1)}\right)\times_{m}\left(U^{(m)}\right)^{T}\right\|_{F}^{2},$
(70)

which is the CP-dictionary-learning problem in [42]. Here $\times_{m}$ denotes
the mode-$m$ product (see [36]) and the outer product is given by

$\operatorname{Out}\left(U^{(1)},\ldots,U^{(m)}\right):=\left[\bigotimes_{k=1}^{m}U^{(k)}[:,1],\bigotimes_{k=1}^{m}U^{(k)}[:,2],\ldots,\bigotimes_{k=1}^{m}U^{(k)}[:,R]\right]\in\mathbb{R}^{I_{1}\times\cdots\times
I_{m}\times R}.$ (71)

The constrained CP-dictionary-learning problem (70) falls under the framework
of block minimization. In fact, let

$\displaystyle\mathbf{A}$
$\displaystyle=\operatorname{Out}\left(U_{n-1}^{(1)},\ldots,U_{n-1}^{(i-1)},U_{n-1}^{(i+1)},\ldots,U_{n-1}^{(m-1)}\right)^{(m)}\in\mathbb{R}^{\left(I_{1}\times\cdots\times
I_{i-1}\times I_{i+1}\times\cdots\times I_{m}\right)\times R}$ (72)
$\displaystyle B$
$\displaystyle=\operatorname{unfold}(\mathbf{A},m)\in\mathbb{R}^{\left(I_{1}\cdots
I_{i-1}I_{i+1}\cdots I_{m}\right)\times R},$ (73)

where unfold(·, $i$ ) denotes the mode- $i$ tensor unfolding (see [36]). Then
the marginal objective function can be rewritten as

$f_{n}^{(i)}(U^{(i)})=\left\|\operatorname{unfold}(X,i)-B\left(U^{(i)}\right)^{T}\right\|_{F}^{2}.$
(74)

Alternating block minimization of (69), i.e. directly minimizing (74) in each
iteration, is known as alternating least squares (ALS), which is studied in
[36, 47, 33, 53] under Euclidean settings, i.e. each $\Theta$ is a convex
subset of $\mathbb{R}^{I_{i}\times R}$. Recently, ALS with diminishing radius
is also studied in [39].

Besides the standard Euclidean setting, we also study the more interesting
setting when some blocks are Riemannian manifolds, e.g. Stiefel manifolds or
rank-$R$ manifolds. In order to solve the corresponding Riemannian CP-
dictionary learning problem (69), instead of directly minimizing
$f_{n}^{(i)}$, we may apply RBMM, which cyclically minimizes a surrogate
function $g_{n}^{(i)}$ in each iteration $n$ given by

$g_{n}^{(i)}(U^{(i)}):=f_{n}^{(i)}(U^{(i)})+\lambda_{n}\|U^{(i)}-U^{(i)}_{n-1}\|_{F}^{2}=\left\|\operatorname{unfold}(X,i)-B\left(U^{(i)}\right)^{T}\right\|_{F}^{2}+\lambda_{n}\|U^{(i)}-U^{(i)}_{n-1}\|_{F}^{2}.$
(75)

As aforementioned at the beginning of Section 5, the Euclidean distance
function is not $g$-smooth with respect to the low-rank manifold. Hence when
at least one $\Theta^{(i)}$ is a low-rank manifold, then in order to ensure
$g$-convexity of the objective function, we will need to take the underlying
manifold to be the Euclidean space, i.e.,
$\mathcal{M}^{(i)}=\mathbb{R}^{I_{i}\times R}$, and allow the constraint set
$\Theta^{(i)}$ to be an embedded submanifold of $\mathbb{R}^{I_{i}\times R}$.
Note that $\Theta^{(i)}$ is not geodesically convex in
$\mathbb{R}^{I_{i}\times R}$ with respect to the Euclidean geometry in
$\mathbb{R}^{I_{i}\times R}$, so we cannot apply our asymptotic complexity
results (Theorems 3.5 and 3.6). However, we can still apply Theorem 3.3 (see
also Remark 3.4). When all the $\Theta^{(i)}$ are either a Stiefel manifold or
the Euclidean space, then we can apply the results in Corollary 3.7 and
therefore have the complexity of $\widetilde{O}(\varepsilon^{-2})$ along with
the asymptotic convergence result.

######  Corollary 5.3 (Complexity of block Euclidean proximal updates for CP-
dictionary learning).

Let $(U_{n}^{(1)},\dots,U_{n}^{(m)})_{n\geq 0}$ be a sequence of factor
matrices obtained by RBMM (Alg. 1) with Euclidean proximal surrogates
$g_{n}^{(i)}$ as in (75). If $\inf_{n}\lambda_{n}>0$, then the sequence of
iterates converges to the set of stationary points of (69). Moreover, if all
the $\Theta^{(i)}$ are either a Stiefel manifold or the Euclidean space, then
the iteration complexity is $\widetilde{O}(\varepsilon^{-2})$.

In the following numerical experiments, we present the advantage of using RBMM
in both the Euclidean setting and the Riemannian setting. We generate
synthetic data $X\in\mathbb{R}^{30\times 20\times 10}$, then apply ALS ([36]),
RBMM with constant $\lambda_{n}=0.1$ and RBMM with diminishing
$\lambda_{n}=0.1\times 0.5^{n}$ as suggested in [47] to find loading matrices
$U^{(1)},U^{(2)},U^{(3)}$ with $R=3$ columns. We run each algorithm 100 times
from the independent random initialization and then compute the averaged
relative error. The typical results of the Euclidean setting, i.e.
$\mathcal{M}^{(i)}=\mathbb{R}^{I_{i}\times R}$ for $i=1,2,3$ are shown in
Figure 6 (A) and (B), where (A) is the case when two RBMM algorithms are
significantly better than ALS. When RBMM reaches a relative error of order
$10^{-14}$, the relative error of ALS is still of order $10^{-3}$; (B) shows
RBMM with constant $\lambda$ is better than ALS while RBMM with diminishing
$\lambda$ is much better than the other two. After 0.25s elapsed time, RBMM
with diminishing $\lambda$ reaches a relative error of order $10^{-14}$, while
that of the other two are still of order $10^{-2}$. After 0.35s elapsed time,
RBMM with constant $\lambda$ reaches error of order $10^{-8}$ while that of
ALS remains $10^{-2}$.

Figure 6 (C) shows the result when $U^{(1)}$ is a point on the Stiefel
manifold $\mathcal{V}^{I_{1}\times R}$ (see Example 4.1), i.e. we let
$\Theta^{(1)}=\mathcal{V}^{I_{1}\times R}$,
$\mathcal{M}^{(2)}=\mathbb{R}^{I_{2}\times R}$ and
$\mathcal{M}^{(3)}=\mathbb{R}^{I_{3}\times R}$. In contrast to the Euclidean
case, all three algorithms demonstrate good performance for random synthetic
data with random initialization. One possible explanation of this phenomenon
is that since $U^{(1)}_{n}\in\mathcal{V}^{I_{1}\times R}$ for all $n$, we have
$(U^{(1)}_{n})^{T}U^{(1)}_{n}=I_{R}$. As a result, the quadratic terms added
in RBMM would reduce to a linear term that does not significantly accelerate
convergence. A similar observation is made when dealing with Stiefel manifolds
in the geodesic subspace tracking problem (see Section 5.1 for details).

Figure 6 (D) shows the convergence results for some synthetic data when
$U^{(1)}$ is a point on fixed-rank matrices manifold $\mathcal{R}_{r}$ with
$r<R$ (see Example 4.2) where $r=5$ and $R=10$, i.e. we let
$\Theta^{(1)}=\mathcal{R}_{r}$, $\mathcal{M}^{(2)}=\mathbb{R}^{I_{2}\times R}$
and $\mathcal{M}^{(3)}=\mathbb{R}^{I_{3}\times R}$. There, RBMM with constant
or diminishing $\lambda$ performs much better than ALS. RBMM with constant and
diminishing $\lambda$ require only $5.16$s and $5.51$s of elapsed time
respectively to reach a relative error of order $10^{-13}$ while that of ALS
takes $9.66$ seconds.

###  5.4. Robust PCA

Principal component analysis (PCA) is a popular technique for reducing the
dimensionality of a data set, while preserving the maximum amount of
information, which is done by seeking the best low-rank approximation of the
high dimensional data set [34]. Mathematically, consider the data matrix
$M\in\mathbb{R}^{m\times n}$, PCA seeks the best $r$-rank approximation of $M$
where $r<\min\\{m,n\\}$, by solving
$\min_{L}\\{\|L-M\|,\operatorname{rank}(L)\leq r\\}$, where $\|X\|$ denotes
the spectral norm. However, it is well known that classical PCA would break
down when entries of $M$ are grossly corrupted. The PCA problem with corrupted
entries is called robust PCA (RPCA) problem [23]. More specifically, it
considers the matrix $M$ of the form $M=L_{0}+S_{0}$ where $L_{0}$ is the
ideal low-rank matrix and $S_{0}$ is a sparse matrix. This problem can be
mathematically formulated as

$\min_{L,S}\,\,\operatorname{rank}(L)+\lambda\|S\|_{0},\quad\text{subject
to}\quad M=L+S,$ (76)

where $\lambda>0$ is a trade-off parameter and $\|S\|_{0}$ is the number of
none zero entries of $S$.

Since the objective function in (76) is the sum of two nonconvex functions, it
is natural to consider a convex relaxation of (76) by relaxing the rank
function by the nuclear norm and the $\ell_{0}$-norm by the $\ell_{1}$-norm.
This gives the following Principal Component Pursuit (PCP) [23, 46]:

$\min_{L,S}\,\,\|L\|_{*}+\lambda\|S\|_{1},\quad\text{subject to}\quad M=L+S$
(77)

where $\|X\|_{*}$ is the nuclear norm of $X$, defined to be the sum of the
singular values. Under certain conditions on the unknown parameters $L_{0}$
and $S_{0}$, it is possible to recover $L_{0}$ and $S_{0}$ by solving the
relaxed problem (77). A more general and realistic setting is that the
observations are noisy, i.e. $M=L+S+Z$ where $Z$ represents the noise term.
This is the stable PCP (SPCP) problem proposed in [63], where the authors
showed the following optimization problem efficiently recovers the true $L$
and $S$ with properly chosen $\mu$,

$\min_{L,S}\,\,\|L\|_{*}+\lambda\|S\|_{1}+\frac{1}{2\mu}\|M-L-S\|_{F}^{2}$
(78)

More studies on solving (78) can also be found in [61, 1].

While (78) and (77) are well-studied formulations of RPCA, they do not always
guarantee a solution with a fixed rank since they promote the rank of $L$ to
be small indirectly through the nuclear norm penalization. Rodriguez and
Wohlberg [55] proposed the following alternative formulation of RPCA by using
a hard low-rank constraint on $L$ instead of the soft nuclear-norm
penalization:

$\min_{L\in\mathcal{R}_{r},\;S\in\mathbb{R}^{m\times
n}}\lambda\|S\|_{1}+\frac{1}{2\mu}\|M-L-S\|_{F}^{2}$ (79)

where $\mathcal{R}_{r}$ is the rank-$r$ matrix manifold (see Example 4.2). In
[55], an alternating minimization algorithm for solving (79) was demonstrated
to be effective through experiments, but no theoretical guarantee of such an
algorithm was provided.

Here, we revisit the alternating minimization algorithm for the Riemannian
formulation of RPCA (79) in [55] through our general framework of RBMM. In
order to apply our framework, we consider a smooth approximation of (79). Note
the 1-norm $\|\cdot\|_{1}$ is in general not smooth. A smooth approximation of
$g(S)=\lambda\|S\|_{1}$ with smoothness parameter $\sigma$ is given by the
following (see e.g. [29])

$g_{\sigma}(S)=\max\left\\{\langle
S,Z\rangle-\frac{\sigma}{2}\|Z\|^{2}_{F}:\|Z\|_{\infty}\leq\lambda\right\\}.$
(80)

It is easy to check the optimal solution of (80) is $Z_{\sigma}(S)$ with

$Z_{\sigma}(S)_{ij}=\min\left\\{\lambda,\max\\{S_{ij}/\sigma,-\lambda\\}\right\\}.$
(81)

Also, the gradient of $g_{\sigma}(S)$ is given by $\nabla
g_{\sigma}(S)=Z_{\sigma}(S)$ and is Lipschitz continuous with parameter
$1/\sigma$.

Now consider the following problem, which is a smooth approximation of the
Riemannian RPCA (79):

$\min_{L\in\mathcal{R}_{r},\;S\in\mathbb{R}^{m\times
n}}F(L,S):=g_{\sigma}(S)+\frac{1}{2\mu}\|M-L-S\|_{F}^{2}.$ (82)

The alternating minimization algorithm for (82) analogous to the one in [55]
for (79) reads as follows. For any $k=1,\cdots,T$, denoting the solution at
iteration $k$ as $L_{k}$ and $S_{k}$, the next iterates are computed by

$\displaystyle L_{k+1}$
$\displaystyle=\operatorname*{arg\,min}_{L\in\mathcal{R}_{r}}F(L,S_{k})=\operatorname*{arg\,min}_{L}\,\,g_{\sigma}(S_{k})+\frac{1}{2\mu}\|M-L-
S_{k}\|_{F}^{2},$ (83) $\displaystyle S_{k+1}$
$\displaystyle=\operatorname*{arg\,min}_{S\in\mathbb{R}^{m\times
n}}F(L_{k+1},S)=\operatorname*{arg\,min}_{S}\,\,g_{\sigma}(S)+\frac{1}{2\mu}\|M-L_{k+1}-S\|_{F}^{2}.$
(84)

The first-order optimality condition of (83) gives $L_{k+1}=U\Sigma_{r}V^{T}$,
where $U\Sigma V^{T}$ is the SVD of $M-S_{k}$ with the singular values ordered
in nonincreasing order and $\Sigma_{r}$ keeps the first $r$ singular values
(see Example 4.2).

Similarly, the first-order optimal solution for (84) is

$(S_{k+1})_{ij}=B_{ij}-\mu\min\left\\{\rho,\max\left\\{-\lambda,\frac{B_{ij}}{\sigma+\mu}\right\\}\right\\}\text{
for }i=1,\ldots,m\text{ and }j=1,\ldots,n$ (85)

where $B=M-L_{k+1}$.

Note that the alternating minimization algorithm in [55] is a special instance
of our RBMM with zero majorization gap. Hence can generalize it by using
proximal regularization as below:

$\displaystyle
L_{k+1}=\operatorname*{arg\,min}_{L\in\mathcal{R}_{r}}G^{(1)}_{k+1}(L):$
$\displaystyle=F(L,S_{k})+\frac{\lambda_{k+1}}{2}\|L-L_{k}\|^{2}_{F}$ (86)
$\displaystyle=g_{\sigma}(S_{k})+\frac{1}{2\mu}\|M-L-
S_{k}\|_{F}^{2}+\frac{\lambda_{k+1}}{2}\|L-L_{k}\|^{2}_{F},$ (87)
$\displaystyle S_{k+1}=\operatorname*{arg\,min}_{S\in\mathbb{R}^{m\times
n}}G^{(2)}_{k+1}(S):$
$\displaystyle=F(L_{k+1},S)+\frac{\lambda_{k+1}}{2}\|S-S_{k}\|^{2}_{F}$ (88)
$\displaystyle=g_{\sigma}(S)+\frac{1}{2\mu}\|M-L_{k+1}-S\|_{F}^{2}+\frac{\lambda_{k+1}}{2}\|S-S_{k}\|^{2}_{F}.$
(89)

The minimization of $G^{(1)}_{k+1}$ and $G^{(2)}_{k+1}$ admits similar closed
form solutions as of (83) and (84), which can be found in Appendix C.4.

We remark again that these iterative updates fall in the framework of RBMM.
Namely, since
$\mathcal{M}^{(1)}=\mathcal{M}^{(2)}=\Theta^{(2)}=\mathbb{R}^{m\times n}$ and
$\Theta^{(1)}=\mathcal{R}_{r}$, the assumptions (A0) and (A1) are satisfied
trivially. Therefore, we deduce the following corollary of Theorem 3.3:

######  Corollary 5.4 (Convergence of BMM for PCP).

Fix a matrix $M\in\mathbb{R}^{m\times n}$. Let
$\boldsymbol{\theta}_{k}:=(L_{k},S_{k})$ be generated by (86) with arbitrary
initialization
$\boldsymbol{\theta}_{0}\in\boldsymbol{\Theta}:=\mathcal{R}_{r}\times\mathbb{R}^{m\times
n}$. Suppose the proximal regularization parameters $\lambda_{k}$ are strictly
positive for all $k\geq 0$. For the smoothed PCP problem (82), the updates in
(86) converge to the set of stationary points.

We remark that in [55], the authors numerically verified the efficiency of the
alternating minimization for solving (79), but no theoretical results on
asymptotic convergence or iteration complexity are established. Corollary 5.4
provides the asymptotic convergence of the algorithm (86) by using the general
RBMM framework.

##  6\. Convergence analysis

In this section, we prove the convergence and complexity of RBMM (Algorithm 1)
stated in Theorems 3.2, 3.3, 3.5, and 3.6. Throughout this section, we let
$(\boldsymbol{\theta}_{n})_{n\geq 1}$ denote the sequence of iterates in
$\boldsymbol{\Theta}\subseteq\prod_{i=1}^{m}\mathcal{M}^{(i)}$ generated by
Algorithm 1. We will also use the notations introduced in Section 2.2. Note in
this section, we use the exponential map as the retraction in the definition
of lifted constraint set (5), i.e.

$\displaystyle T^{*}_{x}\mathcal{M}$ $\displaystyle:=\\{u\in
T_{x}\mathcal{M}\,|\,\textup{$\operatorname{Exp}_{x}(u)=x^{\prime}$ for some
$x^{\prime}\in\boldsymbol{\Theta}$ with $d(x,x^{\prime})\leq r_{0}/2$}\\}.$
(90)

It is worth mentioning that though in the analysis, we will be using an
exponential map, which is computationally expensive, in practice access to it
is not required, as shown in Algorithm 1.

###  6.1. Preliminary analysis

In this section, we present some preliminary results that will be used in the
analysis of RBMM with both $g$-smooth and proximal surrogates.

We first establish basic monotonicity properties of the iterates
$(\boldsymbol{\theta}_{n})_{n\geq 1}$.

######  Proposition 6.1 (Monotonicity of objective and Stability of iterates).

Suppose (A0), except for the $g$-smoothness of $f$. Then the following hold:

(i)

$f(\boldsymbol{\theta}_{n-1})-f(\boldsymbol{\theta}_{n})\geq\sum_{i=1}^{m}\left(g^{(i)}_{n}(\theta^{(i)}_{n})-f^{(i)}_{n}(\theta^{(i)}_{n})-\Delta_{n}\right)\geq-m\Delta_{n}$;

(ii)

$\sum_{n=1}^{\infty}\sum_{i=1}^{m}\left(g^{(i)}_{n}(\theta^{(i)}_{n})-f^{(i)}_{n}(\theta^{(i)}_{n})\right)<f(\boldsymbol{\theta}_{0})-f^{*}+m\sum_{n=1}^{\infty}\Delta_{n}<\infty$.

(iii)

Further assume (A1)(iii) holds. Then the following also holds:

$\displaystyle\sum_{n=1}^{\infty}\sum_{i=1}^{m}\phi\left(d\left(\theta_{n-1}^{(i)},\theta_{n}^{(i)}\right)\right)$
$\displaystyle\leq\sum_{n=1}^{\infty}\sum_{i=1}^{m}\left(g^{(i)}_{n}(\theta^{(i)}_{n})-f^{(i)}_{n}(\theta^{(i)}_{n})\right)<f(\boldsymbol{\theta}_{0})-f^{*}+m\sum_{n=1}^{\infty}\Delta_{n}<\infty.$
(91)

In particular, $d(\theta_{n-1}^{(i)},\theta_{n}^{(i)})=o(1)$ for all
$i=1,\dots,m$.

###### Proof.

Fix $i\in\\{1,\dots,m\\}$. Since $\theta_{n}^{(i\star)}$ is a minimizer of
$g_{n}^{(i)}$ over $\Theta^{(i)}$ (see (A0)(ii)), we get
$g^{(i)}_{n}(\theta^{(i\star)}_{n})\leq
g^{(i)}_{n}(\theta^{(i)}_{n-1})=f^{(i)}_{n}(\theta^{(i)}_{n-1})$, for $n\geq
1$. Hence we deduce

$\displaystyle f^{(i)}_{n}(\theta^{(i)}_{n-1})-f^{(i)}_{n}(\theta^{(i)}_{n})$
$\displaystyle=g^{(i)}_{n}(\theta^{(i)}_{n-1})-g^{(i)}_{n}(\theta^{(i\star)}_{n})+g^{(i)}_{n}(\theta^{(i\star)}_{n})-g^{(i)}_{n}(\theta^{(i)}_{n})+g^{(i)}_{n}(\theta^{(i)}_{n})-f^{(i)}_{n}(\theta^{(i)}_{n})$
(92)
$\displaystyle\geq-\Delta_{n}+g^{(i)}_{n}(\theta^{(i)}_{n})-f^{(i)}_{n}(\theta^{(i)}_{n})\geq-\Delta_{n}.$
(93)

It follows that

$\displaystyle f(\boldsymbol{\theta}_{n-1})-f(\boldsymbol{\theta}_{n})$ (94)
$\displaystyle\qquad=\sum_{i=1}^{m}f(\theta_{n}^{(1)},\dots,\theta_{n}^{(i-1)},\theta_{n-1}^{(i)},\theta_{n-1}^{(i+1)},\dots,\theta_{n-1}^{(m)})-f(\theta_{n}^{(1)},\dots,\theta_{n}^{(i-1)},\theta_{n}^{(i)},\theta_{n-1}^{(i+1)},\dots,\theta_{n-1}^{(m)})$
(95)
$\displaystyle\qquad=\sum_{i=1}^{m}f_{n}^{(i)}(\theta^{(i)}_{n-1})-f_{n}^{(i)}(\theta^{(i)}_{n})$
(96)
$\displaystyle\qquad\geq\sum_{i=1}^{m}g^{(i)}_{n}(\theta^{(i)}_{n})-f^{(i)}_{n}(\theta^{(i)}_{n})-\Delta_{n}$
(97) $\displaystyle\qquad\geq-m\Delta_{n}.\quad$ (98)

This shows (i).

Next, to show (ii), adding up the above inequality for $n\geq 1$,

$\displaystyle\sum_{n=1}^{\infty}\sum_{i=1}^{m}\left(g^{(i)}_{n}(\theta^{(i)}_{n})-f^{(i)}_{n}(\theta^{(i)}_{n})\right)$
$\displaystyle\leq\left(\sum_{n=1}^{\infty}f(\boldsymbol{\theta}_{n-1})-f(\boldsymbol{\theta}_{n})\right)+m\sum_{n=1}^{\infty}\Delta_{n}\leq
f(\boldsymbol{\theta}_{0})-f^{*}+m\sum_{n=1}^{\infty}\Delta_{n}<\infty.$ (99)

To show (iii), note that by (A1)(iii),

$g^{(i)}_{n}(\theta^{(i)}_{n})-f^{(i)}_{n}(\theta^{(i)}_{n})\geq\phi\left(d\left(\theta_{n-1}^{(i)},\theta_{n}^{(i)}\right)\right).$
(100)

Together with (ii) we get (iii). ∎

The following is an immediate consequence of Proposition 6.1 and sub-level
compactness.

######  Proposition 6.2 (Boundedness of iterates).

Under (A0) except for the $g$-smoothness of $f$, the set
$\\{\boldsymbol{\theta}_{n}:n\geq 1\\}$ is bounded.

###### Proof.

Let $m\sum_{n=1}^{\infty}\Delta_{n}=T<\infty$, by Prop. 6.1 we immediately
have $f(\boldsymbol{\theta}_{n})\leq
f(\boldsymbol{\theta}_{0})+m\sum_{k=1}^{n}\Delta_{k}<f(\boldsymbol{\theta}_{0})+T$.
Consider
$K=\\{\boldsymbol{\theta}\in\boldsymbol{\Theta}:f(\boldsymbol{\theta})\leq
f(\boldsymbol{\theta}_{0})+T\\}$, by (A0)(i), $K$ is compact. Hence the set
$\\{\boldsymbol{\theta}_{n}:n\geq 1\\}$ is bounded. ∎

The following two propositions to be introduced will only be used in deriving
the rate of convergence results, and they assume the geodesic convexity of the
constraint sets.

If $\mathcal{M}$ is a Riemannian manifold and if $x,y\in\mathcal{M}$ such that
$d(x,y)<\operatorname{r_{\textup{inj}}}(x)$, then the tangent vector
$\eta_{x}(y):=\operatorname{Exp}_{x}(y)\in T_{x}\mathcal{M}$ is uniquely
defined. Now we introduce the following notations:

$\displaystyle\eta_{n+1}^{(i)}(n):=\eta_{\theta_{n+1}^{(i)}}(\theta_{n}^{(i)}),\quad\eta_{n}^{(i)}(n+1):=\eta_{\theta_{n}^{(i)}}(\theta_{n+1}^{(i)}).$
(101)

Note that the above tangent vectors are well-defined for all $n$ sufficiently
large under (A0) and (A1)(i)-(iii). Indeed, by Proposition 6.1 we have
$d(\boldsymbol{\theta}_{n},\boldsymbol{\theta}_{n+1})=o(1)$. Also, the
injectivity radius is uniformly lower bounded by some constant $r_{0}>0$ (see
(A1)(ii)). Hence

$\displaystyle
d(\boldsymbol{\theta}_{n},\boldsymbol{\theta}_{n+1})<r_{0}\quad\textup{for all
$n\geq N_{0}$ for some $N_{0}\in\mathbb{N}$},$ (102)

so the tangent vectors in (101) are well-defined for $n\geq N_{0}$.

First, we show that the first-order variation of the objective function along
the trajectory of iterates $\boldsymbol{\theta}_{n}$ is summable.

######  Proposition 6.3 (Finite first-order variation).

Suppose (A0), (A1)(i)-(ii), and (A1)(iii) with $\phi(x)=x^{2}$ hold. Further
assume that the constraint set $\Theta^{(i)}$ is strongly convex in
$\mathcal{M}^{(i)}$ for $i=1,\dots,m$. Then

$\sum_{n=N_{0}}^{\infty}\sum_{i=1}^{m}\left\lvert\left\langle\operatorname{grad}f_{n}^{(i)}(\theta_{n+1}^{(i)}),\eta_{n+1}^{(i)}(n)\right\rangle\right\rvert<\frac{L_{f}+2}{2}\left(f(\boldsymbol{\theta}_{0})-f^{*}\right)+3m\sum_{n=0}^{\infty}\Delta_{n+1}<\infty.$
(103)

###### Proof.

By Lemma B.1,

$\left\lvert
f_{n}^{(i)}(\theta_{n}^{(i)})-f_{n}^{(i)}(\theta_{n+1}^{(i)})-\left\langle\operatorname{grad}f_{n}^{(i)}(\theta_{n+1}^{(i)}),\eta^{(i)}_{n+1}(n)\right\rangle\right\rvert\leq\frac{L_{f}}{2}d^{2}(\theta_{n}^{(i)},\theta_{n+1}^{(i)}).$
(104)

Then by Proposition 6.1,

$\displaystyle\left\lvert\left\langle\operatorname{grad}f_{n}^{(i)}(\theta_{n+1}^{(i)}),\eta^{(i)}_{n+1}(n)\right\rangle\right\rvert$
$\displaystyle\leq\frac{L_{f}}{2}d^{2}(\theta_{n}^{(i)},\theta_{n+1}^{(i)})+\left\lvert
f_{n}^{(i)}(\theta_{n}^{(i)})-f_{n}^{(i)}(\theta_{n+1}^{(i)})\right\rvert$
(105)
$\displaystyle\leq\frac{L_{f}}{2}d^{2}(\theta_{n}^{(i)},\theta_{n+1}^{(i)})+f_{n}^{(i)}(\theta_{n}^{(i)})-f_{n}^{(i)}(\theta_{n+1}^{(i)})+2\Delta_{n+1}.$
(106)

Summing up for all $i=1,\dots,m$ and then for $n\geq 1$, we get

$\displaystyle\sum_{n=N_{0}}^{\infty}\sum_{i=1}^{m}\left\lvert\left\langle\operatorname{grad}f_{n}^{(i)}(\theta_{n+1}^{(i)}),\eta^{(i)}_{n+1}(n)\right\rangle\right\rvert$
(107)
$\displaystyle\qquad\leq\frac{L_{f}}{2}\sum_{n=N_{0}}^{\infty}d^{2}(\boldsymbol{\theta}_{n},\boldsymbol{\theta}_{n+1})+\sum_{n=N_{0}}^{\infty}\left(f(\boldsymbol{\theta}_{n})-f(\boldsymbol{\theta}_{n+1})\right)+2m\sum_{n=N_{0}}^{\infty}\Delta_{n+1}$
(108)
$\displaystyle\qquad\leq\frac{L_{f}}{2}\sum_{n=0}^{\infty}\phi(d(\boldsymbol{\theta}_{n},\boldsymbol{\theta}_{n+1}))+f(\boldsymbol{\theta}_{0})-f^{*}+2m\sum_{n=0}^{\infty}\Delta_{n+1}$
(109)
$\displaystyle\qquad\leq\frac{L_{f}+2}{2}\left(f(\boldsymbol{\theta}_{0})-f^{*}\right)+3m\sum_{n=0}^{\infty}\Delta_{n+1}$
(110)

where the last inequality comes from Proposition 6.1.

∎

Next, we relate the first-order approximation of the objective difference
$f(\boldsymbol{\theta}_{n+1})-f(\boldsymbol{\theta}_{n})$ with that of the
worst-case difference $f(\boldsymbol{\theta})-f(\boldsymbol{\theta}_{n})$.
This will be used crucially in the proof of Theorem 3.5 and 3.6.

######  Proposition 6.4 (Asymptotic first-order optimality).

Suppose (A0) and (A1)(i)-(iii) hold. Further assume that the constraint set
$\Theta^{(i)}$ is strongly convex in $\mathcal{M}^{(i)}$ for $i=1,\dots,m$.
Fix a sequence $(b_{n})_{n\geq 1}$ of positive reals with
$b_{n}\leq\hat{r}:=\min\\{r_{0},1\\}$ for all $n$. Then the following hold for
all $n\geq N_{0}$:

(i)

($g$-smooth surrogates) Suppose (A1)(i-gs) holds. Then there exists constant
$c_{1},c_{2}>0$ independent of $\boldsymbol{\theta}_{0}$ such that

$\displaystyle\sum_{i=1}^{m}\left\langle\operatorname{grad}f_{n+1}^{(i)}(\theta_{n}^{(i)}),\eta^{(i)}_{n}(n+1)\right\rangle$
(111)
$\displaystyle\qquad\leq\hat{r}^{-1}b_{n+1}\sum_{i=1}^{m}\inf_{\eta^{(i)}\in
T^{*}_{\theta_{n}^{(i)}},\|\eta^{(i)}\|\leq
1}\left\langle\operatorname{grad}g_{n+1}^{(i)}(\theta_{n}^{(i)})+\operatorname{grad}_{i}f(\boldsymbol{\theta}_{n})-\operatorname{grad}f_{n+1}^{(i)}(\theta^{(i)}_{n}),\eta^{(i)}\right\rangle$
(112)
$\displaystyle\qquad\qquad+c_{1}b_{n+1}^{2}+c_{2}d^{2}(\boldsymbol{\theta}_{n},\boldsymbol{\theta}_{n+1})+m\Delta_{n}.$
(113)

(ii)

(Riemannian proximal surrogates) Suppose (A1)(i-rp) holds. Then there exists
constant $c_{1},c_{2}>0$ independent of $\boldsymbol{\theta}_{0}$ such that

$\displaystyle\sum_{i=1}^{m}\left\langle\operatorname{grad}f_{n+1}^{(i)}(\theta_{n}^{(i)}),\eta^{(i)}_{n}(n+1)\right\rangle$
$\displaystyle\leq\hat{r}^{-1}b_{n+1}\sum_{i=1}^{m}\inf_{\eta^{(i)}\in
T^{*}_{\theta_{n}^{(i)}},\|\eta^{(i)}\|\leq
1}\left\langle\operatorname{grad}f_{n+1}^{(i)}(\theta_{n}^{(i)}),\eta^{(i)}\right\rangle$
(114)
$\displaystyle\qquad+c_{1}b_{n+1}^{2}+c_{2}d^{2}(\boldsymbol{\theta}_{n},\boldsymbol{\theta}_{n+1})+m\Delta_{n}.$
(115)

###### Proof.

We first show (ii). Let $\hat{r}=\min\\{r_{0},1\\}$, where $r_{0}$ is the
lower bound of injectivity radius (see (A1) (ii)). Fix arbitrary
$\boldsymbol{\theta}=[\theta^{(1)},\cdots,\theta^{(m)}]\in\boldsymbol{\Theta}$
such that

$\displaystyle d(\theta^{(i)},\theta_{n}^{(i)})\leq
b_{n+1}\leq\hat{r}\leq\operatorname{r_{\textup{inj}}}(\theta^{(i)})$ (116)

for all $i$. For each $i$, let
$\gamma^{(i)}:[0,1]\rightarrow\mathcal{M}^{(i)}$ be the unique distance-
minimizing geodesic from $\theta_{n}^{(i)}$ to $\theta^{(i)}$. Denote
$\eta^{(i)}:=(\gamma^{(i)})^{\prime}(0)$. Then $\eta^{(i)}\in
T^{*}_{\theta_{n}^{(i)}}$. Denote $h_{n}^{(i)}:=g_{n}^{(i)}-f_{n}^{(i)}$.

Fix $i=1,\dots,m$. We will drop the superscripts $(i)$ indicating block $i$ in
the proof until we mention otherwise. As in (A0)(ii), let
$\theta_{n+1}^{\star}$ be an exact minimizer of $g_{n+1}$ over $\Theta$. Then

$\displaystyle\langle\operatorname{grad}f_{n+1}(\theta_{n}),\eta_{n}(n+1)\rangle-\frac{L_{f}}{2}d^{2}(\theta_{n},\theta_{n+1})$
(117)
$\displaystyle\qquad\overset{(a)}{\leq}f_{n+1}(\theta_{n+1})-f_{n+1}(\theta_{n})$
(118)
$\displaystyle\qquad\overset{(b)}{\leq}g_{n+1}(\theta_{n+1})-g_{n+1}(\theta_{n})$
(119)
$\displaystyle\qquad\overset{(c)}{\leq}g_{n+1}(\theta^{\star}_{n+1})+\Delta_{n}-g_{n+1}(\theta_{n})$
(120)
$\displaystyle\qquad\overset{(d)}{\leq}f_{n+1}(\gamma(1))-f_{n+1}(\theta_{n})+\Delta_{n}+h_{n+1}(\gamma(1))$
(121)
$\displaystyle\qquad\overset{(e)}{\leq}\langle\operatorname{grad}f_{n+1}(\theta_{n}),\eta\rangle+\frac{L_{f}}{2}d^{2}(\theta_{n},\theta)+\Delta_{n}+h_{n+1}(\gamma(1)).$
(122)

Here, $(a)$ follows from geodesic $L_{f}$-smoothness of $f$ (see (A0)(i)),
$(b)$ follows from definition of majorizing surrogates (Def. 2.1), $(c)$
follows from the definition of $\theta_{n}^{\star}$ and the optimality gap
$\Delta_{n}$ (see (A0)(ii)), $(d)$ follows from the fact that
$g_{n+1}(\theta_{n})=f_{n+1}(\theta_{n})$ and
$g_{n+1}(\theta_{n+1}^{\star})\leq g_{n+1}(\gamma(t))$ for all $t\in[0,1]$,
$(e)$ uses geodesic $L_{f}$-smoothness of $f$. Also,

$\displaystyle
h_{n+1}(\gamma(1))=\frac{\lambda_{n}}{2}d^{2}(\theta,\theta_{n})\leq\frac{\lambda_{n}}{2}b_{n+1}^{2}.$
(123)

This gives

$\displaystyle\left\langle\operatorname{grad}f_{n+1}(\theta_{n}),\,\eta_{n}(n+1)\right\rangle$
$\displaystyle\leq\left\langle\operatorname{grad}f_{n+1}(\theta_{n}),\,\eta\right\rangle+c_{1}b_{n+1}^{2}+c_{2}d^{2}(\theta_{n},\theta_{n+1})+\Delta_{n},$
(124)

where $c_{1}:=\frac{L_{f}+\lambda_{n}}{2}$ and $c_{2}=\frac{L_{f}}{2}$.

The above holds for all $\theta\in\Theta$ with $d(\theta,\theta_{n})\leq
b_{n+1}$. Hence

$\displaystyle\left\langle\operatorname{grad}f_{n+1}(\theta_{n}),\,\eta_{n}(n+1)\right\rangle$
$\displaystyle\leq\inf_{\|\eta\|\leq b_{n+1},\eta\in
T^{*}_{\theta_{n}}}\left\langle\operatorname{grad}f_{n+1}(\theta_{n}),\,\eta\right\rangle+c_{1}b_{n+1}^{2}+c_{2}d^{2}(\theta_{n},\theta_{n+1})+\Delta_{n}.$
(125)

Now observe by g-convexity of $\Theta$, if $\alpha\eta\in T^{*}_{\theta_{n}}$
for $\alpha\geq 1$, then $\eta\in T^{*}_{\theta_{n}}$. Therefore we get

$\displaystyle\sup_{\|\eta\|\leq\hat{r},\eta\in
T^{*}_{\theta_{n}}}\left\langle-\operatorname{grad}f_{n+1}(\theta_{n}),\,\eta\right\rangle$
$\displaystyle=\sup_{\|u\|\leq b_{n+1},\hat{r}b^{-1}_{n+1}u\in
T^{*}_{\theta_{n}}}\left\langle-\operatorname{grad}f_{n+1}(\theta_{n}),\,\hat{r}b^{-1}_{n+1}u\right\rangle$
(126) $\displaystyle\leq\hat{r}b^{-1}_{n+1}\sup_{\|u\|\leq b_{n+1},u\in
T^{*}_{\theta_{n}}}\left\langle-\operatorname{grad}f_{n+1}(\theta_{n}),\,u\right\rangle$
(127)

The above together with (125) show (ii).

Lastly, we note that (i) be shown by a similar argument. The main difference
is that we do not have an explicit expression for the majorization gap
$h_{n+1}$ as in (123) but instead, we have $g$-smoothness of the surrogates.
Hence we can proceed by using the $g$-smoothness of $g_{n+1}^{(i)}$ in line
(120) instead of decomposing $g_{n+1}^{(i)}=f_{n+1}^{(i)}+h_{n+1}^{(i)}$ and
using the $g$-smoothness of $f_{n+1}^{(i)}$. The rest of the analysis is
identical. ∎

###  6.2. RBMM with Riemannian proximal surrogates

In this section, We analyze our RBMM algorithm with Riemannian proximal
surrogates. Under (A1)(i-rp), Algorithm 1 is the Riemannian analog of the
Euclidean block majorization-minimization with Euclidean proximal regularizer.

Before we present the proof of our general results in Theorem 3.3 and 3.5, we
will first delve into the special case of block Riemannian proximal methods on
Hadamard manifolds, as briefly mentioned in Section 4.3. Examining this
special case will be helpful for understanding RBMM with Riemannian proximal
surrogates, particularly in two key aspects. First, the Riemannian proximal
surrogate $g_{n}^{(i)}$ is geodesically strongly convex on Hadamard manifolds,
which makes the subproblem for the minimization step in (LABEL:eq:RBMM_MmMm)
easy to solve. Second, $d^{2}(\theta,\theta^{(i)}_{n-1})$ is not $g$-smooth
even on Hadamard manifolds, highlighting the necessity of analyzing the case
of $g$-smooth surrogates and Riemannian/Euclidean proximal surrogates
separately.

We first recall some useful facts on Hadamard manifolds. First, recall that on
Hadamard manifolds, the exponential map is a global diffeomorphism, and
therefore its injectivity radius is $+\infty$ [56, Theorem 4.1, p.221].
Second, $d^{2}(\cdot,p):\mathcal{M}\to\mathbb{R}$ is geodesically $2$-strongly
convex for fixed $p\in\mathcal{M}$ where $\mathcal{M}$ is a Hadamard manifold
([14]). We give the formal definition of geodesically strongly convex below.

######  Definition 6.5 (Geodesic strong convexity).

Suppose $\mathcal{M}$ is a Hadamard manifold. A function
$f:\mathcal{M}\to\mathbb{R}$ is geodesically $\alpha$-strongly convex if it
satisfies

$f(y)-f(x)-\left\langle\operatorname{grad}f(x),\operatorname{Exp}_{x}^{-1}(y)\right\rangle_{x}\geq\frac{\alpha}{2}d^{2}(x,y)$
(128)

for all $x$, $y\in\mathcal{M}$, where $d(x,y)$ is the geodesic distance
between $x$ and $y$ on $\mathcal{M}$.

We can also show by definition, that the Riemannian proximal surrogate
$g_{n}^{(i)}(\theta)$ is geodesically strongly convex with appropriate
Riemannian proximal regularization parameter. This is formally stated in the
following proposition.

######  Proposition 6.6 (Geodesical strong convexity of surrogates).

Let $\mathcal{M}$ be a Hadamard manifold and $f:\mathcal{M}\to\mathbb{R}$ is
geodesically $L_{f}$-smooth (see Definition 3.1). Fix $p\in\mathcal{M}$, let
$g:\mathcal{M}\to\mathbb{R}$

$g(x)=f(x)+\frac{c}{2}d^{2}(x,p)$ (129)

where $c>L_{f}$. Then $g$ is geodesically strongly convex with parameter
$c-L_{f}$.

###### Proof.

First by $g$-smoothness of $f$ using Lemma B.1 and the geodesical 2-strong
convexity of $d^{2}(\cdot,p)$,

$\displaystyle
f(y)-f(x)-\left\langle\operatorname{grad}f(x),\operatorname{Exp}_{x}^{-1}(y)\right\rangle_{x}$
$\displaystyle\geq-\frac{L_{f}}{2}d^{2}(x,y)$ (130) $\displaystyle
d^{2}(y,p)-d^{2}(x,p)-\left\langle\operatorname{grad}d^{2}(x,p),\operatorname{Exp}_{x}^{-1}(y)\right\rangle_{x}$
$\displaystyle\geq d^{2}(x,y)$ (131)

for any $x$, $y\in\mathcal{M}.$ Finally multiplying the second inequality by
$\frac{c}{2}$ and adding up with the first inequality finish the proof. ∎

Hence, a major benefit of using Riemannian proximal surrogates on Hadamard
manifolds is that the sub-problems in RBMM are geodesically strongly convex so
can be efficiently solved by using standard convex Riemannian optimization
methods (see, e.g., [58, 64, 41]).

Another benefit from the geodesic strong convexity of surrogates is related to
the inexact computation, which is stated in the following proposition. Namely,
the square of the geodesic distance between the inexact and the exact
minimizer of the proximal surrogate (20) is upper bounded by the optimality
gap (18), when the surrogate is geodesically strongly convex. Therefore, (19)
in (A0)(ii) is automatically satisfied on Hadamard manifolds.

######  Proposition 6.7 (Optimality gap for iterates).

For each $n\geq 1$ and $i\in\\{1,\ldots,m\\}$, let $\theta_{n}^{(i\star)}$ be
the exact minimizer of the geodesically
$\left(\lambda_{n}-L_{f}\right)$-strongly convex function $\theta\mapsto
g_{n}^{(i)}(\theta)$ in (20) over the strongly convex set $\Theta^{(i)}$. Then

$\frac{\lambda_{n}-L_{f}}{2}d^{2}\left(\theta_{n}^{(i\star)},\theta_{n}^{(i)}\right)\leq\Delta_{n}.$

###### Proof.

Note the optimality condition of $\theta_{n}^{(i\star)}$ over $\Theta^{(i)}$
is

$\langle
g_{n}^{(i)}(\theta_{n}^{(i\star)}),\operatorname{Exp}^{-1}_{\theta_{n}^{(i\star)}}(\theta)\rangle\geq
0,\quad\forall\theta\in\Theta^{(i)}$ (132)

then we have

$\displaystyle\frac{\lambda_{n}-L_{f}}{2}d^{2}\left(\theta_{n}^{(i\star)},\theta_{n}^{(i)}\right)$
$\displaystyle\leq
g_{n}^{(i)}(\theta_{n}^{(i)})-g_{n}^{(i)}(\theta_{n}^{(i\star)})-\langle
g_{n}^{(i)}(\theta_{n}^{(i\star)}),\operatorname{Exp}^{-1}_{\theta_{n}^{(i\star)}}(\theta_{n}^{(i)})\rangle$
(133) $\displaystyle\leq
g_{n}^{(i)}(\theta_{n}^{(i)})-g_{n}^{(i)}(\theta_{n}^{(i\star)})$ (134)
$\displaystyle\leq\Delta_{n}.$ (135)

∎

Next, we give an illustration of $d^{2}(\cdot,p)$ is not $g$-smooth in
general. The following proposition gives the expression of the Riemannian
gradient of geodesic distance.

######  Proposition 6.8 (Riemannian gradient of geodesic distance).

Let $\mathcal{M}$ be a complete Riemannian manifold, $p\in\mathcal{M}$ with
$r_{\operatorname{inj}}(p)=r$. Consider the function
$h:\mathcal{M}\rightarrow\mathbb{R}$ where $h(x)=d_{\mathcal{M}}^{2}(x,p)$,
where $p\in\mathcal{M}$ is fixed. If $d(x,p)<r$, then
$\operatorname{grad}(h)=-2\operatorname{Exp}_{x}^{-1}(p)$ as a vector in
$T_{x}\mathcal{M}$.

###### Proof.

See e.g. [56, Proposition 4.8, p. 108]. ∎

With Prop. 6.8, we can now explain why the surrogate given in (A1)(i-rp) is
not $g$-smooth in general. Recall Def. 3.1, in order $\frac{1}{2}d^{2}(x,p)$
in Prop. 6.8 to be $g$-smooth, we need the following to hold,

$\frac{1}{2}\left\lVert\operatorname{grad}d^{2}(y,p)-\Gamma_{x\rightarrow
y}(\operatorname{grad}d^{2}(x,p))\right\rVert=\left\lVert-\operatorname{Exp}_{y}^{-1}(p)+\Gamma_{x\rightarrow
y}(\operatorname{Exp}_{x}^{-1}(p))\right\rVert\leq Ld(x,y)$ (136)

In Euclidean space, parallel transport $\Gamma_{x\rightarrow y}$ is identical,
and (136) would degenerate to equality with $L=1$. Hence
$\frac{1}{2}d^{2}(x,p)$ is $g$-smooth in Euclidean space, see Fig. 7(a) for an
illustrative example. However, on a general Riemannian manifold, even on the
Hadamard manifold, this is no longer true in general, see Fig. 7(b) for a
counterexample of hyperbolic space (see Example 4.4). Some more examples of
$S^{1}$ are shown in Fig. 8.

Figure 7.  Examples on $g$-smoothness of $d^{2}(x,p)$. Panel (a) is an example
in Euclidean space, the length of the dashed black arrow represents the LHS of
(136) and the length of the black segment is $d(x,y)$; Panel (b) is a
counterexample in hyperbolic space.  Figure 8.  Examples on $g$-smoothness of
$d^{2}(x,p)$ on $S^{1}$. Panel (a) is the case when (136) does not hold; Panel
(b) (c) is the case when (136) becomes an equality with $L=2$.

Now we come back to establish the proof of Theorems 3.3 and 3.5. Let
$\boldsymbol{\theta}_{\infty}=[\theta^{(1)}_{\infty},\cdots,\theta^{(m)}_{\infty}]$
be a limit point of $(\boldsymbol{\theta}_{n})_{n}$, for simplicity of
notations, replace $(\boldsymbol{\theta}_{n})_{n}$ by the convergent
subsequence. We introduce the following notations of tangent vectors and sets,
which will be used throughout the remaining proofs. For $n\geq 1$, let

$\displaystyle\Theta^{(i)}_{n}=\\{\theta\in\Theta^{(i)}:d(\theta,\theta^{(i)}_{n})\leq
r_{0}/2\\}$ (137)

where $r_{0}$ is the uniform lower bound of injectivity radius (see (A1)(ii)).
Similarly, denote

$\Theta^{(i)}_{\infty}=\\{\theta\in\Theta^{(i)}:d(\theta,\theta^{(i)}_{\infty})\leq
r_{0}/2\\}.$ (138)

For any $\theta^{(i)}\in\Theta^{(i)}_{\infty}$, let

$\displaystyle\eta_{n}^{(i)}:=\eta_{\theta_{n}^{(i)}}(\theta^{(i)}),\quad\eta_{n}^{(i\star)}:=\eta_{\theta_{n}^{(i\star)}}(\theta^{(i)}),\;\textup{and}\;\;\eta_{\infty}^{(i)}:=\eta_{\theta_{\infty}^{(i)}}(\theta^{(i)}).$
(139)

Note that the above tangent vectors are well-defined for all $n$ sufficiently
large under (A0) and (A1)(i)-(ii). In fact, since $\boldsymbol{\theta}_{n}$
converges to $\boldsymbol{\theta}_{\infty}$,

$\displaystyle
d(\theta^{(i)}_{n},\theta^{(i)}_{\infty})<r_{0}/4\quad\textup{for all $n\geq
N_{1}$ for some $N_{1}\in\mathbb{N}$}.$ (140)

Also, by (A0),

$\displaystyle
d(\theta^{(i)}_{n},\theta^{(i\star)}_{n})<r_{0}/4\quad\textup{for all $n\geq
N_{2}$ for some $N_{2}\in\mathbb{N}$}.$ (141)

Therefore by triangle inequality, the tangent vectors in (139) are well
defined for $n\geq\max\\{N_{1},N_{2}\\}$. The following proposition shows the
asymptotic first-order optimality of inexact solutions, which is the key
proposition to prove Theorem 3.3.

######  Proposition 6.9 (Asymptotic first order optimality of inexact
solutions).

Assume (A1)(i-rp), (ii), (A0). Let
$\boldsymbol{\theta}_{\infty}=[\theta^{(1)}_{\infty},\cdots,\theta^{(m)}_{\infty}]$
be a limit point of $(\boldsymbol{\theta}_{n})_{n}$, for simplicity of
notations, replace $(\boldsymbol{\theta}_{n})_{n}$ by the convergent