Inexact Adaptive Cubic Regularization Algorithms on Riemannian Manifolds and Application

We consider the large-scale separable unconstrained optimization problem on general Riemannian manifolds:

where $\mathcal{M}$ is a Riemannian manifold, $n\gg 1$, and each $f_{i}:\mathcal{M}\rightarrow\mathbb{R}$ is continuously differentiable ($i:=1,2,\dots n$). Such problems frequently appear in machine learning and scientific computing, where each $f_{i}$ is a loss (or misfit) function corresponding to $i-$th observation (or measurement); see, e.g., [1, 2, 3]. In such “large-scale” settings, since the evaluations of the gradient or the Hessian of $f$ can be computationally expensive, some inexact techniques are used for approximating the first and the second derivatives, which particularly includes the random sampling technique. As a result, many popular first-order and second-order stochastic algorithms are proposed to solve problem (1). The stochastic gradient descent (SGD) algorithm is one of the first-order stochastic/inexact algorithms, which uses only the gradient information. The idea of the SGD comes from [4] and one can find some variants in [5, 6, 27, 10, 11] and the references therein (the Riemannian versions of the SGD can be found in [10, 11]). The second-order stochastic/inexact algorithms always use the inexact Hessian information, which includes the popular sub-sampling Newton method, sub-sampling/inexact trust region algorithm, and the sub-sampling/inexact cubic regularization algorithm; see e.g, [28, 29, 30, 31, 32, 17] (the Riemannian version of the inexact trust region algorithm is in [17]). Particularly, we note the inexact trust-region algorithm and the inexact adaptive cubic regularization algorithms introduced in [32] for solving problem (1). Both algorithms employ the true gradient and the inexact Hessian (by sub-sampling) of $f$ and solve the corresponding sub-problems approximately every iteration. These two algorithms are respectively shown to own iteration complexity $\mathcal{O}(\max\{\varepsilon_{g}^{-2}\varepsilon_{H}^{-1},\varepsilon_{H}^{-3}\})$ and $\mathcal{O}(\max\{\varepsilon_{g}^{-2},\varepsilon_{H}^{-3}\})$ of achieving the $(\varepsilon_{g},\varepsilon_{H})$-optimality (see definition of the $(\varepsilon_{g},\varepsilon_{H})$-optimality in Definition 2.2).

It is known that the Riemannian geometry framework has some advantages in many applications, such as translating some nonconvex (constrained) optimization problems in Euclidean space into convex (unconstrained) ones over a Riemannian manifold; see, e.g., [12, 34, 35, 36]. As a result, some classical numerical methods for solving optimization problems on the Euclidean space, such as Newton’s method, BFGS algorithm, trust region method, gradient algorithm, subgradient algorithm, etc., have been extended to the Riemannian manifold setting; see, e.g., [16, 12, 13, 14, 15, 34, 33, 35, 36]. As we have noted above, for solving large-scale problem (1), the stochastic gradient descent algorithm and its variants have been extended from Euclidean spaces to Riemannian manifolds in the literature; see, e.g., [10, 11]. Recently, the Riemannian inexact trust region algorithm is studied in [17], which uses the inexact gradient and the inexact Hessian every iteration, and in addition the inexact solution of the trust region sub-problem. Noting that the algorithm employs the inexact gradient instead of the true gradient, it extends the corresponding one in [32] to more general cases even in Euclidean settings. Similarly, for achieving the $(\varepsilon_{g},\varepsilon_{H})$-optimality, the iteration complexity of the algorithm is also proven to be $\mathcal{O}(\max\{\varepsilon_{g}^{-2}\varepsilon_{H}^{-1},\varepsilon_{H}^{-3}\})$. At the same time, some numerical results are provided which illustrate the algorithm performs significantly better than state-of-the-art deterministic and stochastic algorithms in some applications.

Inspired by the prior works (in particular, the works in [32] and [17]), we propose the Riemannian inexact adaptive cubic regularization algorithm for solving problem (1) on general Riemannian manifolds, which employs the inexact gradient and the inexact Hessian every iteration. The algorithm is proven to own the similar iteration complexity $\mathcal{O}(\max\{\varepsilon_{g}^{-2},\varepsilon_{H}^{-3}\})$ for obtaining the $(\varepsilon_{g},\varepsilon_{H})$-optimality under certain conditions on accuracies about the inexact gradient and the inexact Hessian. Particularly, iteration complexities of the Riemannian (deterministic) adaptive cubic regularization algorithm and the Riemannian inexact adaptive cubic regularization algorithm under the true gradient are established. As an application, the proposed algorithms are applied to solve the joint diagonalization problem on the Stiefel manifold. Numerical results indicate that inexact algorithms are more efficient than the deterministic algorithm at the same accuracy. Meanwhile, the inexact Riemannian adaptive cubic regularization algorithm outperforms the inexact Riemannian trust-region algorithm (in [17, Algorithm 1]).

The paper is organized as follows. As usual, some basic notions and notation on Riemannian manifolds, together with some related properties about the sub-sampling, are introduced in the next section. The Riemannian inexact adaptive cubic regularization algorithm and its iteration complexity are presented in section 3, and the application to joint diagonalization on the Stiefel manifold is shown in the last section.

Notation and terminologies used in the present paper are standard; the readers are referred to some textbooks for more details; see, e.g., [12, 18, 19, 22].

Let $\mathcal{M}$ be a connected and complete $n$-dimensional Riemannian manifold. Let $x\in\mathcal{M}$, and let ${T}_{x}\mathcal{M}$ stand for the tangent space at $x$ to $\mathcal{M}$. $T\mathcal{M}:=\cup_{x\in\mathcal{M}}T_{x}\mathcal{M}$ is called the tangent bundle on $\mathcal{M}$. We denote by $\langle,\rangle_{x}$ the scalar product on $T_{x}\mathcal{M}$ with the associated norm $\|\cdot\|_{x}$, where the subscript $x$ is sometimes omitted. For $y\in\mathcal{M}$, let $\gamma:[0,1]\rightarrow\mathcal{M}$ be a piecewise smooth curve joining $x$ to $y$. Then, the arc-length of $\gamma$ is defined by $l(\gamma):=\int_{0}^{1}\|{\gamma}^{\prime}(t)\|dt$, while the Riemannian distance from $x$ to $y$ is defined by ${\rm d}(x,y):=\inf_{\gamma}l(\gamma)$, where the infimum is taken over all piecewise smooth curves $\gamma:[0,1]\rightarrow\mathcal{M}$ joining $x$ to $y$.

We use $\nabla$ to denote the Levi-Civita connection on $\mathcal{M}$. A vector field $V$ is said to be parallel along $\gamma$ if $\nabla_{{\gamma}^{\prime}}V=0$. In particular, for a smooth curve $\gamma$, if ${\gamma}^{\prime}$ is parallel along itself, then $\gamma$ is called a geodesic; thus, a smooth curve $\gamma$ is a geodesic if and only if $\nabla_{{\gamma}^{\prime}}{{\gamma}^{\prime}}=0$. A geodesic $\gamma:[0,1]\rightarrow M$ joining $x$ to $y$ is minimal if its arc-length equals its Riemannian distance between $x$ and $y$. By the Hopf-Rinow theorem [18], $(\mathcal{M},{\rm d})$ is a complete metric space, and there is at least one minimal geodesic joining $x$ to $y$ for any points $x$ and $y$.

Let $f:\mathcal{M}\rightarrow{{\mathbb{R}}}$ be a twice continuously differentiable real-valued function defined on $\mathcal{M}$. The Remannian gradient and Hessian of $f$ at $x\in\mathcal{M}$ are denoted by ${\rm grad}f(x)$ and ${\rm Hess}f(x)$, which are respectively defined as

On manifolds $\mathcal{M}$, we use the retraction defined as below to move in the direction of a tangent vector, which can be found in [12, Definition 4.1.1 and Proposition 5.5.5].

Let $J:=\{1,2,\cdots,n\}$. We consider the large-scale minimization (1):

where $n\gg 1$, and for each $i\in J$, $f_{i}:\mathcal{M}\rightarrow\mathbb{R}$ is a twice continuously differentiable real-valued function. Note that $f$ may be non-convex, some authors refer to find a proximate $(\varepsilon_{g},\varepsilon_{H})$-optimality in practice, which is defined by [32, Definition 1] (see also [17, Definition 2.1]). As usual, $\lambda_{\min}(A)$ stands the minimum eigenvalue of matrix $A$.

In practice, we adopt the sub-sampled inexact gradient and Hessian. To proceed, let $S_{g},S_{H}\subset J$ be the sample collections with or without replacement from $J$, and their cardinalities are denoted as $|S_{g}|$ and $|S_{H}|$, respectively. The sub-sampled inexact gradient and Hessian are defined as

respectively. The following lemma provides sufficient sample sizes of $S_{g}$ and $S_{H}$ to guarantee that the inexact $G$ and $H$ approximate ${\rm grad}f$ and ${\rm Hess}f$ in a probabilistic way, which is taken form [17, Theorem 4.1]. Here we set

In this section, we propose the following inexact Riemannian adaptive cubic regularization algorithm for solving problem (1), which employs the inexact Hessian and gradient at each iteration.

Let $\{x_{k}\}$ be a sequence generated by Algorithm 3.1. We make the following blanket assumptions throughout the paper:

As stated in Step 4, we solve the sub-problem (3) approximately every iteration. The most popular conditions for the approximate solution in the literature are the Cauchy and Eigenpoint conditions; see, e.g., [23, 24, 32]. To ensure the convergence of the algorithm, we make the following similar assumptions on the approximate solutions $\{\eta_{k}\}$:

Now, let $k\in\mathbb{N}$. Noting that the sub-problem (3) of Algorithm 3.1 is actually posed on Euclidean spaces, the follow properties of $\eta_{k}^{C}$ and $\eta_{k}^{E}$ are valid (see [32, Lemmas 6,7]).

Lemma 3.2 below estimates the sufficient decrease in the objective function.

In the following lemma, we estimate the upper bound for the cubic regularization parameters $\{\sigma_{k}\}$ before the algorithm terminates. Here all iterations satisfying (12) are said to be successful and to be unsuccessful otherwise.

Without loss generality, we assume that $\sigma_{0}\leq 2\gamma L_{H}$. Moreover, from now on, let $\mathcal{N}_{succ}$ and $\mathcal{N}_{fail}$ stand the set of all the successful and unseccessful iterations of Algorithm 3.1, respectively. The number of successful (unsuccessful) iterations is denoted by $|\mathcal{N}_{succ}|$ $(|\mathcal{N}_{fail}|)$. The lemma below provides estimation of the total number of successful iterations before the algorithm terminates.

We can now prove the main theorem, which indicates the iteration complexity of Algorithm 3.1.

As corollaries of Theorem 3.6, we consider the (deterministic) Riemannian adaptive cubic regularization algorithm and the inexact Riemannian adaptive cubic regularization algorithm under the true gradient. The following corollary shows the iteration complexity of the (deterministic) Riemannian adaptive cubic regularization algorithm, which is immediate from Theorem 3.6 (noting that assumption (A3) is satisfied naturally).

Following the same line of proof of Theorem 3.6, we have the following corollary about the iteration complexity of the inexact Riemannian adaptive cubic regularization algorithm under the true gradient.