Oracle complexities of augmented Lagrangian methods for nonsmooth manifold optimization

When the objective function is geodesically convex over a Riemannian manifold, the subgradient-type methods are studied in [4, 15, 14]. An asymptotic convergence result is first established in [14], while an iteration complexity of $\mathcal{O}(\epsilon^{-2})$ is established in [4, 15]. For the general nonsmooth problem on the Riemannian manifold, the Riemannian subgradient may not exist. By sampling points within a neighborhood of the current iteration point at which the objective function $f$ is differentiable, the Riemannian gradient sampling methods are developed in [18, 19].

There are many works addressing nonsmooth problems in the form of (1.1). When $\mathcal{A}$ in (1.1) is the identity operator and the proximal operator of $h$ can be calculated analytically, Riemannian proximal gradient methods [10, 22, 23] are proposed. These methods achieve an outer iteration complexity of $\mathcal{O}(\epsilon^{-2})$ for attaining an $\epsilon$-stationary solution. It should be emphasized that each iteration in these algorithms involves solving a subproblem that lacks an explicit solution, and they utilize the semismooth Newton method to solve it. For the general $\mathcal{A}$, Peng et al. [37] develop the Riemannian smoothing gradient method by utilizing a smoothing technique from the Moreau envelope. They show that the algorithm achieves an iteration complexity of $\mathcal{O}(\epsilon^{-3})$. Moreover, Beck and Rosset [3] propose a dynamic smoothing gradient descent on manifold and obtain an iteration complexity of $\mathcal{O}(\epsilon^{-3})$. Very recently, based on a similar smoothing technique, a Riemannian alternating direction method of multipliers (RADMM) is proposed in [29] with an iteration complexity result of $\mathcal{O}(\epsilon^{-4})$ for driving a Karush–Kuhn–Tucker (KKT) residual based stationarity. This complexity result is a bit worse than $\mathcal{O}(\epsilon^{-3})$ [37] because there is a lack of adaptive adjustment of the smoothing parameter. There is also a line of research based on the augmented Lagrangian function. Deng and Peng [12] proposed a manifold inexact augmented Lagrangian method for solving problem (1.1). Zhou et al. [46] consider a more general problem and propose a globalized semismooth Newton method to solve its subproblem efficiently. In these two papers, they only give the asymptotic convergence result, while the iteration complexity is missing.

In the case when the manifold $\mathcal{M}$ is specified by equality constraints $c(x)=0$, e.g., the Stiefel manifold. Problem (1.1) can be regarded as a constrained optimization problem with a nonsmooth and nonconvex objective function. Given that $\mathcal{M}$ is often nonconvex, we only list the related works in case of that the constraint functions are nonconvex. Papers [32, 38] propose and study the iteration complexity of augmented Lagrangian methods for solving nonlinearly constrained nonconvex composite optimization problems. The iteration complexity results they achieve for an $\epsilon$-stationary point are both $\tilde{\mathcal{O}}(\epsilon^{-3})$. More specifically, [38] uses the accelerated gradient method of [16] to obtain the approximate stationary point. On the other hand, the authors in [32] obtain such approximate stationary point by applying an inner accelerated proximal method as in [9, 26], whose generated subproblems are convex. It is worth mentioning that both of these papers make a strong assumption about how the feasibility of an iterate is related to its stationarity. Lin et al. [34] propose an inexact proximal-point penalty method by solving a sequence of penalty subproblems. Under a non-singularity condition, they show a complexity result of $\tilde{\mathcal{O}}(\epsilon^{-3})$.

Finally, we point out that there exist several works for nonsmooth expectation problems on the Riemannian manifold, i.e., when $f(x)=\mathbb{E}[f(x,\xi)]$ in problem (1.1). In particular, Li et al. [30] focus on a class of weakly convex problems on the Stiefel manifold. They present a Riemannian stochastic subgradient method and show that it has an iteration complexity of $\mathcal{O}(\epsilon^{-4})$ for driving a natural stationarity measure below $\epsilon$. Peng et al. [37] propose a Riemannian stochastic smoothing method and give an iteration complexity of $\mathcal{O}(\epsilon^{-5})$ for driving a natural stationary point. Li et al. [33] propose a stochastic inexact augmented Lagrangian method to solve problems involving a nonconvex composite objective and nonconvex smooth functional constraints. Under certain regularity conditions, they establish an oracle complexity result of $\mathcal{O}(\epsilon^{-5})$ with a stochastic first-order oracle. Wang et al. [40] propose two Riemannian stochastic proximal gradient methods for minimizing a nonsmooth function over the Stiefel manifold. They show that the proposed algorithm finds an $\epsilon$-stationary point with $\mathcal{O}(\epsilon^{-3})$ first-order oracle complexity. It should be emphasized that each oracle call involves a subroutine whose iteration complexity remains unclear.

In Table 1, we summarize our complexity results and several existing methods to produce an $\epsilon$-stationary point. The definitions of first-order oracle in determined and stochastic settings are provided in Definition 3.12 and 4.9, respectively. The notation $\tilde{\mathcal{O}}$ denotes the dependence of a $\log$ factor. Here, we only consider the scenario where both the objective function and the nonlinear constraint are nonconvex. For the iteration complexity of ALM in convex case, please refer to [2, 35, 43, 42, 27], etc. We do not add the proximal type algorithms in [10, 40, 22, 23] since those algorithms involve using the semismooth Newton method to solve the subproblem, where the overall oracle complexity is not clear. It can easily be shown that our algorithms achieve better oracle complexity results both in determined and stochastic settings. In particular, in the stochastic setting, our algorithm establishes an oracle complexity of $\tilde{\mathcal{O}}(\epsilon^{-3.5})$, which is better than the best-known $\mathcal{O}(\epsilon^{-4})$ result.

Let $\left<\cdot,\cdot\right>$ and $\|\cdot\|$ be the Euclidean product and induced Euclidean norm. Given a matrix $A$, we use $\|A\|_{F}$ to denote the Frobenius norm, $\|A\|_{1}:=\sum_{ij}|A_{ij}|$ to denote the $\ell_{1}$ norm. For a vector $x$, we use $\|x\|_{2}$ and $\|x\|_{1}$ to denote its Euclidean norm and $\ell_{1}$ norm, respectively. The indicator function of a set $\mathcal{C}$, denoted by $\delta_{\mathcal{C}}$, is set to 0 on $\mathcal{C}$ and $+\infty$ otherwise. The distance from $x$ to $\mathcal{C}$ is denoted by $\mathrm{dist}(x,\mathcal{C}):=\min_{y\in\mathcal{C}}\|x-y\|$. We use $\nabla f(x)$ and $\operatorname{grad}f(x)$ to denote the Euclidean gradient and Riemannian gradient of $f$, respectively.

An $n$-dimensional smooth manifold $\mathcal{M}$ is an $n$-dimensional topological manifold equipped with a smooth structure, where each point has a neighborhood that is diffeomorphism to the $n$-dimensional Euclidean space. For all $x\in\mathcal{M}$, there exists a chart $(U,\psi)$ such that $U$ is an open set and $\psi$ is a diffeomorphism between $U$ and an open set $\psi(U)$ in the Euclidean space. A tangent vector $\eta_{x}$ to $\mathcal{M}$ at $x$ is defined as tangents of parametrized curves $\gamma$ on $\mathcal{M}$ such that $\gamma(0)=x$ and

where $\wp_{x}\mathcal{M}$ is the set of all real-valued functions $f$ defined in a neighborhood of $x$ in $\mathcal{M}$. Then, the tangent space $T_{x}\mathcal{M}$ of a manifold $\mathcal{M}$ at $x$ is defined as the set of all tangent vectors at point $x$. The manifold $\mathcal{M}$ is called a Riemannian manifold if it is equipped with an inner product on the tangent space $T_{x}\mathcal{M}$ at each $x\in\mathcal{M}$. If $\mathcal{M}$ is a Riemannian submanifold of an Euclidean space $\mathcal{E}$, the inner product is defined as the Euclidean inner product: $\left<\eta_{x},\xi_{x}\right>=\mathrm{tr}(\eta_{x}^{\top}\xi_{x})$. The Riemannian gradient $\operatorname{grad}f(x)\in T_{x}\mathcal{M}$ is the unique tangent vector satisfying

If $\mathcal{M}$ is a compact Riemannian manifold embedded in an Euclidean space, we have that $\operatorname{grad}f(x)=\mathcal{P}_{T_{x}\mathcal{M}}(\nabla f(x))$, where $\nabla f(x)$ is the Euclidean gradient, $\mathcal{P}_{T_{x}\mathcal{M}}$ is the projection operator onto the tangent space $T_{x}\mathcal{M}$. The retraction operator is one of the most important ingredients for manifold optimization, which turns an element of $T_{x}\mathcal{M}$ into a point in $\mathcal{M}$.

We have the following Lipschitz-type inequalities on the retraction on the compact submanifold.

Another basic ingredient for manifold optimization is the vector transport.

We first introduce the concept of the Moreau envelope.

The following result demonstrates that the Moreau envelope of a convex function is not only continuously differentiable but also possesses a bounded gradient norm.

Next, we provide the definition of retraction smoothness. This concept plays a crucial role in the convergence analysis of the algorithms proposed in the subsequent sections.

Similar to the Lipschitz continuity of gradient in Euclidean space, we give the following extended definition of the Lipschitz continuity of the Riemannian gradient.

The following lemma shows that the Lipschitz continuity of $\operatorname{grad}f(x)$ can be deduced by the Lipschitz continuity of $\nabla f(x)$.

We say $x\in\mathcal{M}$ is an $\epsilon$-stationary point of problem (1.1) if there exists $\zeta\in\partial h(\mathcal{A}x)$ such that $\left\|\mathcal{P}_{T_{x}\mathcal{M}}\left(\nabla f(x)+\mathcal{A}^{*}\zeta\right)\right\|\leq\epsilon$. As in [45], there is no finite time algorithm that can guarantee $\epsilon$-stationarity in the nonconvex nonsmooth setting. Motivated by the notion of $(\delta,\epsilon)$-stationarity introduced in [45], we introduce a generalized $\epsilon$-stationarity for problem (1.1) based on the Lagrangian function.

Since problem (1.1) contains the nonsmooth term, we introduce the auxiliary variable $y=\mathcal{A}x$, which results in the following equivalent split problem:

The Lagrangian function $l:\mathcal{M}\times\mathbb{R}^{m}\times\mathbb{R}^{m}\rightarrow\mathbb{R}$ of (3.1) is given by:

where $z$ is the corresponding Lagrangian multiplier. Then the KKT condition of (1.1) is given as follows [12]: A point $x\in\mathcal{M}$ satisfies the KKT condition if there exists $y,z\in\mathbb{R}^{m}$ such that

Based on the KKT condition, we give the definition of $\epsilon$-stationary point for (1.1):

Before constructing the augmented Lagrangian function of (1.1), we first focus on the equivalent split problem (3.1). In particular, the augmented Lagrangian function of (3.1) is constructed as follows:

where $\sigma$ is the penalty parameter. It is shown that the augmented Lagrangian method involves the variable $x$ and the auxiliary variables $y$. However, for the original problem (1.1), which only features the variable $x$, the corresponding augmented Lagrangian function should solely include $x$ as well as the associated multipliers.

An intuitive way of introducing the augmented Lagrangian function associated with (1.1) is to eliminate the auxiliary variables in (3.5), i.e.,

Note that the optimal solution $y^{*}$ of the above minimization problem has the following closed-form solution:

Therefore, one can obtain the explicit form of $\mathcal{L}_{\sigma}(x,z)$:

where $M_{h}^{1/\sigma}$ is the Moreau envelope given in Definition 2.4. It is important to emphasize that the augmented Lagrangian function (3.6) only retains $x$ as the primal variable. The Lagrangian multiplier $z$ corresponds to the nonsmooth function $h$. When comparing $\mathcal{L}_{\sigma}$ with $\mathbb{L}_{\sigma}$, it becomes evident that $\mathcal{L}_{\sigma}$ involves fewer variables. We will also demonstrate that the augmented Lagrangian function (3.6) possesses favorable properties, such as being continuously differentiable. This characteristic proves to be advantageous when it comes to designing and analyzing algorithms [13, 31].

As shown in [12], a manifold augmented Lagrangian method based on the augmented Lagrangian function (3.5) can be written as

By the definition of $\mathcal{L}_{\sigma}$ and (3.7), we know that (3.9) can be rewritten as

We call (3.10) the augmented Lagrangian method based on (3.6) for (1.1). The authors in [12] show that an inexact version of (3.10) converges to a stationary point. However, the iteration complexity result is unknown. Motivated by [38], we design a novel manifold inexact augmented Lagrangian framework, as outlined in Algorithm 1 to efficiently find an $\epsilon$-KKT point of (1.1). The key distinction from (3.10) lies in the update strategy for the Lagrangian multiplier $z^{k}$ and the careful selection of $\sigma_{k}$. This strategic choice will enable us to bound the norm of $z^{k}$ and ensures that $\frac{\|z^{k}\|}{\sigma_{k}}$ converges to zero.

In Algorithm 1, we introduce the dual stepsize $\beta_{k}$, which is used in the update of the dual variable $z^{k+1}$ (also called the Lagrangian multiplier). The update rule (3.13) of $\beta_{k+1}$ is to make sure that the dual variable $z^{k+1}$ is bounded when $k\rightarrow\infty$. The detailed derivation is given in the proof of Theorem 3.8. Moreover, for the convenience of checking the termination criteria and subsequent convergence analysis, we introduce an auxiliary sequence denoted as $y^{k}$. The main difference between our algorithm and the one presented in [12] lies in our algorithm’s reliance on a more refined selection of step sizes and improved control over the accuracy of $x$-subproblem solutions. Another crucial difference is the stopping criteria for the subproblem. In particular, we provide two options based on either a given residual or a given number of iterations, further enhancing the flexibility of the algorithm These modifications allow us to derive results pertaining to the iterative complexity of the algorithm.

The function $\mathcal{L}_{\sigma}(x,z)$ with respect to $x$ exhibits favorable properties, such as continuous differentiability, enabling a broader range of choices in selecting subproblem-solving algorithms. Let us denote $\psi_{k}(x):=\mathcal{L}_{\sigma_{k}}(x,z^{k})$. By the property of the proximal mapping and the Moreau envelope in Proposition 2.5, one can readily conclude that $\psi_{k}(x)$ is continuously differentiable and

The main computational cost of Algorithm 1 lies in solving the subproblem (3.11). In the next subsection, we will give a subroutine to solve (3.11).

In this subsection, we give a subroutine to find each $x^{k+1}$ in Algorithm 1. We first present a general algorithm for smooth problems on Riemannian manifolds, and then try to apply the algorithm to solve the subproblem (3.11) in Algorithm 1 by checking the related condition. Let us focus on the following smooth problem on the Riemannian manifold:

where $\mathcal{M}$ is a Riemannian submanifold embedded in Euclidean space. Denote $\psi_{\min}:=\inf_{x\in\mathcal{M}}\psi(x)$. We assume that $\psi$ is retr-smooth with $L$ and $\mathcal{R}$, i.e.,

We adopt a Riemannian gradient descent method (RGD) [1, 6] to solve problem (3.16). The detailed algorithm is given in Algorithm 2. The following lemma gives the iteration complexity result of Algorithm 2.

To apply Algorithm 2 for computing $x^{k+1}$ in Algorithm 1, it is crucial to demonstrate that $\psi_{k}(x)$ is retr-smooth, meaning it satisfies (3.17). The following lemma provides the essential insight.

With Lemmas 3.2 and 3.4, we can apply Algorithm 2 to find each $x^{k+1}$. The lemma below gives the inner iteration complexity for the $k$-th outer iteration of Algorithm 1.

We have shown the oracle complexity for solving the subproblems of ManIAL. Now, let us first investigate the outer iteration complexity and conclude this subsection with total oracle complexity.