A Riemannian Alternating Descent Ascent Algorithmic Framework for Nonconvex-Linear Minimax Problems on Riemannian Manifolds

Meng Xu and Ya-Feng Liu are supported in part by the National Natural Science Foundation of China (NSFC) under Grant 12371314 and Grant 12021001. Bo Jiang is supported by the National Natural Science Foundation of China (NSFC) under Grant 12522116 and Grant 12371314. Anthony Man-Cho So is supported in part by the Hong Kong Research Grants Council (RGC) General Research Fund (GRF) project CUHK 14204823.

In this paper, we consider a class of nonconvex-linear minimax problems on Riemannian manifolds, which find wide applications in machine learning and signal processing. For solving this class of problems, we develop a flexible Riemannian alternating descent ascent (RADA) algorithmic framework. Within this framework, we propose two easy-to-implement yet efficient algorithms that alternately perform one or multiple projected/Riemannian gradient descent steps and a proximal gradient ascent step at each iteration. We show that the proposed RADA algorithmic framework can find both an ε-Riemannian-game-stationary point and an ε-Riemannian-optimization-stationary point within 𝒪⁡(ε−3) iterations, achieving the best-known iteration complexity. We also reveal intriguing similarities and differences between the algorithms developed within our proposed framework and existing algorithms, thus providing important insights into the improved efficiency of the former. Lastly, we present numerical results on sparse principal component analysis (PCA), fair PCA, and sparse spectral clustering to demonstrate the superior performance of the proposed algorithms.

1 Introduction

In this paper, we consider Riemannian nonconvex-linear (NC-L) minimax problems of the form

minx∈ℳmaxy∈ℰ2{F(x,y):=f(x)+⟨𝒜(x),y⟩−h(y)}, (1)

where ℳ is a Riemannian manifold embedded in a finite-dimensional Euclidean space ℰ1, ℰ2 is another finite-dimensional Euclidean space, f:ℰ1→ℝ is a continuously differentiable function, 𝒜:ℰ1→ℰ2 is a smooth mapping, and h:ℰ2→(−∞,+∞] is a proper closed convex function with a compact domain and a tractable proximal mapping. Many applications in machine learning and signal processing give rise to instances of problem (1). We briefly introduce three such applications below and refer the interested reader to [11, 23, 42, 73, 27, 40, 50] for more applications.

1.1 Motivating Applications

In the following, we use A=[a1,a2,…,aN]∈ℝd×N to denote a data matrix, where each of the N columns corresponds to a data sample with d attributes.

Sparse principal component analysis (SPCA). Each principal component of A obtained through the classic PCA is a linear combination of d attributes, making it difficult to interpret the derived principal components as new features [84]. To achieve a good balance between dimension reduction and interpretability, SPCA seeks principal components with very few nonzero components. This motivates the following formulation of SPCA [33]:

minX∈𝒮⁡(d,r)⁡{−⟨A​A⊤,X​X⊤⟩+μ​‖X‖1}. (2)

Here, 𝒮⁡(d,r)={X∈ℝd×r∣X⊤​X=Ir} is the Stiefel manifold with Ir being the r-by-r identity matrix, μ>0 is the weighting parameter, and ‖X‖1=∑i,j|Xi​j| is the ℓ1-norm of the matrix X. By Fenchel duality, the SPCA problem (2) can be equivalently reformulated as

minX∈𝒮⁡(d,r)⁡maxY∈𝒴​{−⟨A​A⊤,X​X⊤⟩+⟨Y,X⟩}, (3)

where 𝒴={Y∈ℝd×r∣‖Y‖∞≤μ} with ‖Y‖∞=maxi,j⁡|Yi​j|. Problem (3) is an instance of problem (1) with ℳ being the Stiefel manifold, 𝒜⁡(X)=X, and h being the indicator function of the compact set 𝒴.

Fair principal component analysis (FPCA). Suppose that the N data samples belong to m groups according to certain clustering, and each group i corresponds to a d×ni submatrix Ai with ∑i=1mni=N. Recall that classic PCA aims to find a subspace with dimension r<d to minimize the total reconstruction error, or equivalently, to maximize the total variance. As such, certain groups may suffer a higher reconstruction error than the others [62]. To reduce such disparity among different groups, FPCA minimizes the maximum reconstruction error among the m groups, which is equivalent to maximizing the minimum variance among the m groups. Mathematically, the FPCA problem can be formulated as [62, 79, 74, 64]

minX∈𝒮⁡(d,r)⁡maxi∈{1, 2,…,m}−⟨Ai​Ai⊤,X​X⊤⟩,

which is equivalent to

minX∈𝒮⁡(d,r)⁡maxy∈Δm−∑i=1myi​⟨Ai​Ai⊤,X​X⊤⟩. (4)

Here, Δm={y∈ℝm|∑i=1myi=1,yi≥0,i=1, 2,…,m} is the standard simplex in ℝm. Note that problem (4) is an instance of problem (1), where ℳ is the Stiefel manifold, 𝒜:ℝd×r→ℝm is given by 𝒜​(X)i=−⟨Ai​Ai⊤,X​X⊤⟩ for i=1,2,…,m, and h is the indicator function of Δm.

Sparse spectral clustering (SSC). This task aims to divide N data samples into m groups, each of which consists of similar data points. Spectral clustering constructs a symmetric affinity matrix W=[Wi​j]N×N, where Wi​j≥0 measures the pairwise similarity between two samples ai and aj. To promote the sparsity and interpretability of spectral clustering, the works [52, 51, 70] proposed SSC, which entails solving the following problem:

minQ∈𝒢⁡(N,m)⁡{⟨L,Q⟩+μ​‖Q‖1}. (5)

Here, 𝒢⁡(N,m)={X​X⊤∈ℝN×N∣X∈𝒮⁡(N,m)} is the Grassmann manifold [6], which can be regarded as a submanifold in the Euclidean space [63]; L=IN−S−1/2WS−1/2 is the normalized Laplacian matrix with S1/2 being the diagonal matrix with diagonal elements s1,s2,…,sN and si=∑jWi​j; and μ>0 is the weighting parameter. By Fenchel duality, problem (5) can be equivalently reformulated as

minQ∈𝒢⁡(N,m)⁡maxY∈𝒴​{⟨L,Q⟩+⟨Y,Q⟩}. (6)

Again, problem (6) is an instance of problem (1) with ℳ being the Grassmann manifold, 𝒜⁡(Q)=Q, and h being the indicator function of the compact set 𝒴:={Y∈ℝN×N∣‖Y‖∞≤μ}.

1.2 Related Works

Minimax problems in the Euclidean space—i.e., the setting where ℳ⊆ℰ1 is an arbitrary closed convex set in problem (1)—have been the focus of many recent studies. For general smooth nonconvex-concave (NC-C) minimax problems, many nested-loop algorithms have been developed and analyzed [56, 68, 45, 60, 78, 36, 57]. Nevertheless, single-loop algorithms have gained growing interest due to their simplicity. A classic single-loop algorithm is the gradient descent ascent (GDA) method, which performs a (projected) gradient descent step on x and a (projected) gradient ascent step on y at each iteration. For smooth nonconvex-strongly concave (NC-SC) minimax problems, GDA can find an ε-stationary point within 𝒪⁡(ε−2) iterations [46]. However, when F⁡(x,⋅) is not strongly concave for some x, GDA may encounter oscillations even for bilinear games. Many improved GDA variants have been proposed to overcome this issue [32, 53, 81, 77, 76, 41, 72, 75, 46]. One powerful approach is to design an appropriate regularized objective function. Such an approach was used by Xu et al. [76] and Zhang et al. [81] to develop the alternating gradient projection algorithm and the smoothed GDA algorithm, respectively. Both algorithms can find an ε-stationary point of a general smooth NC-C minimax problem within 𝒪⁡(ε−4) iterations. Later, Li et al. [41] extended the techniques in [81] and developed the smoothed proximal linear descent ascent algorithm. The algorithm enjoys, among other things, an iteration complexity of 𝒪⁡(ε−4) for finding an ε-stationary point of certain nonsmooth NC-C minimax problem. For the special case of NC-L minimax problems, Pan et al. [58], Shen et al. [65], and He et al. [20] used the regularization approach to develop algorithms that can reach ε-stationarity within 𝒪⁡(ε−3) iterations, which is the best complexity known to date.

Recently, there has been growing interest in minimax problems on Riemannian manifolds, but existing works on the topic are few. These works can be divided into two groups based on the nature of the constraints: (i) The constraints on x and y are manifolds [82, 34, 19, 18] and (ii) the constraint on x is a manifold while the constraint on y is a closed convex set [26, 74]. Among the works in the first group, Zhang et al. [82] generalized the classic Sion’s minimax theorem to geodesic metric spaces and proposed a Riemannian corrected extragradient (RCEG) algorithm for solving geodesically convex-geodesically concave (GC-GC) minimax problems. Jordan et al. [34] analyzed both the last-iterate and average-iterate convergence of RCEG and a Riemannian gradient descent ascent (RGDA) algorithm in various GC-GC settings. For problems without geodesic convexity, Han et al. [19] extended the Hamiltonian gradient method to the Riemannian setting and established its global linear convergence under the assumption that the Riemannian Hamiltonian of F satisfies the Riemannian Polyak-Łojasiewicz (PŁ) condition. In addition, Han et al. [18] proposed several second-order methods that are proven to asymptotically converge to local minimax points of nonconvex-nonconcave minimax problems under strong assumptions. As for the works in the second group, Huang and Gao [26] proposed a RGDA algorithm different from the one in [34] and established its iteration complexity of 𝒪⁡(ε−2) for finding an ε-stationary point of an NC-SC minimax problem. Later, the work [74] presented an alternating Riemannian/projected gradient descent ascent (ARPGDA) algorithm with an iteration complexity of 𝒪⁡(ε−3) for finding an ε-stationary point of the Riemannian NC-L minimax problem (1) when h is the indicator function of a convex compact set.

From the above discussion, we observe that existing works on minimax problems over Riemannian manifolds do not address the settings of problem (1) that we are interested in. Indeed, neither the geodesic convexity assumption in [82, 34] nor the strong concavity assumption in [26] apply to the motivating applications introduced in Section 1.1. Additionally, the Riemannian PŁ condition of the Hamiltonian function in [19] is difficult to verify and is not known to hold for these applications. On the other hand, the second-order methods in [18] cannot be used to tackle problem (1) since the Hessian of F⁡(x,⋅) is singular. Given the above background, we are motivated to develop simple yet efficient single-loop first-order methods with strong theoretical guarantees for solving the Riemannian NC-L minimax problem (1).

One possible approach is based on the following reformulation of (1) as a Riemannian nonsmooth composite problem:

minx∈ℳ⁡{f⁡(x)+h∗​(𝒜⁡(x))}. (7)

Here, h∗ is the conjugate function of h. Indeed, since h is assumed to have a compact domain, we know from [5, Theorem 4.23] that h∗ is Lipschitz continuous. Then, we can apply a host of algorithms to tackle problem (7), including Riemannian subgradient-type methods [7, 22, 21, 42, 25], Riemannian proximal gradient-type methods [11, 27, 28, 70, 12, 49], Riemannian smoothing-type algorithms [4, 59, 80], and splitting-type methods [39, 38, 14, 40, 83, 13, 47]. In particular, when 𝒜 is the identity map ℐ and h∗ is a real-valued weakly convex function, Li et al. [42] proposed a family of Riemannian subgradient (RSG) methods and showed that all of them have an iteration complexity of 𝒪⁡(ε−4) for getting certain stationarity measure of problem (7) below ε. When 𝒜=ℐ, Chen et al. [11] proposed the manifold proximal gradient (ManPG) method with an outer iteration complexity of 𝒪⁡(ε−2) for finding an ε-stationary point of problem (7). Here, the outer iteration complexity refers to the number of subproblems that need to be solved in the nested-loop algorithms. In particular, the complexity of solving each subproblem is not counted in the outer iteration complexity. When 𝒜 is a linear mapping, Beck and Rosset [4] proposed a dynamic smoothing gradient descent method (DSGM) with an iteration complexity of 𝒪⁡(ε−3) for finding an ε-stationary point of problem (7). Using the splitting technique, Li et al. [40] proposed a Riemannian ADMM (RADMM) with an iteration complexity of 𝒪⁡(ε−4) and Deng et al. [13] proposed a manifold inexact augmented Lagrangian (ManIAL) method with an iteration complexity of 𝒪⁡(ε−3) for finding an ε-stationary point of problem (7). When 𝒜 is a nonlinear mapping, Wang et al. [70] proposed a manifold proximal linear (ManPL) method and established its outer iteration complexity of 𝒪⁡(ε−2) for finding an ε-stationary point of problem (7).

In this paper, we propose an alternative approach, which is based on exploiting the minimax structure of problem (1) and can also be used to solve problem (7). As seen from Table 1, which summarizes the applicability and complexity of different methods, our approach is able to address more general problem settings and has the best complexity known to date among single-loop algorithms for solving problem (1). Moreover, as we shall see in Section 5, our approach reveals intriguing similarities and differences between our proposed algorithms for problem (1) and existing algorithms for problem (7).

Table 1: A summary of methods with iteration complexity results for solving problem (1) or (7), where “⋆” denotes the outer iteration complexity. The results for HiBSA (developed for general Euclidean NC-C minimax problems) and AGP (developed for Euclidean NC-L minimax problems) are stated for the Euclidean NC-L minimax setting.

Algorithm Mapping 𝒜 Function h∗ Complexity Type

RSG [42] ℐ real-valued weakly convex 𝒪⁡(ε−4) single-loop

ManPG [11] ℐ convex, Lipschitz continuous 𝒪​(ε−2)⋆ nested-loop

ManIAL [13] linear convex, Lipschitz continuous 𝒪⁡(ε−3) nested-loop

DSGM [4] linear convex, Lipschitz continuous 𝒪⁡(ε−3) single-loop

RADMM [40] linear convex, Lipschitz continuous 𝒪⁡(ε−4) single-loop

ManPL [70] nonlinear convex, Lipschitz continuous 𝒪​(ε−2)⋆ nested-loop

ARPGDA [74] nonlinear support function of a convex compact set 𝒪⁡(ε−3) single-loop

HiBSA [53] nonlinear support function of a convex compact set 𝒪⁡(ε−4) single-loop

AGP [58] nonlinear support function of a convex compact set 𝒪⁡(ε−3) single-loop

RADA-PGD/RGD [this paper] nonlinear convex, Lipschitz continuous 𝒪⁡(ε−3) single-loop

1.3 Our Contributions

We now summarize the main contributions of this paper as follows.

We propose a flexible Riemannian alternating descent ascent (RADA) algorithmic framework for solving problem (1). At the heart of our framework is the observation that the value function associated with the maximization of certain regularized version of the objective function F with respect to y is smooth and has an easily computable gradient; see (11) and (18). At each iteration, the proposed RADA algorithmic framework finds an approximate minimizer of the said value function to update x and performs one proximal gradient ascent step to update y. The value function contains information about the best response of y, namely, the maximization of the objective function F⁡(x,y) with respect to y for any given x. Therefore, minimizing the proposed value function can be regarded as approximately solving a surrogate of problem (1). It should be noted that our approach to updating x is different from simply finding an (approximate) minimizer of some surrogate of F⁡(⋅,y) for a fixed y, as is done in, e.g., [76, 58, 53, 20], and is more advantageous in practice. Within the RADA algorithmic framework, we provide two customized simple yet efficient single-loop first-order algorithms, namely, RADA-PGD and RADA-RGD, which perform a fixed number of projected gradient descent steps and Riemannian gradient descent steps to update x at each iteration, respectively. Interestingly, even when adapted to the setting where ℳ is simply a closed convex set in a Euclidean space, our proposed RADA algorithmic framework, along with the two customized single-loop algorithms, is different from existing methods for solving Euclidean NC-L minimax problems (such as those in [76, 58, 53, 20]) and is, to the best of our knowledge, new.

We prove that our proposed RADA algorithmic framework can find an ε-Riemannian-game-stationary (ε-RGS) point and an ε-Riemannian-optimization-stationary (ε-ROS) point of problem (1) within 𝒪⁡(ε−3) iterations. This iteration complexity result matches the best-known complexity result for general NC-L minimax problems where the variable x is constrained on a nonempty compact convex set instead of a Riemannian manifold [58, 20]. It also matches the best-known complexity result for the Riemannian nonsmooth composite problem (7) with a linear mapping 𝒜 [4, 13]. Furthermore, while one typically needs extra computation to convert between an ε-game stationary point and an ε-optimization stationary point in the Euclidean setting (see, e.g., [46, 77]), it is worth noting that using our RADA algorithmic framework, the complexity of finding an ε-RGS point is of the same order as that of finding an ε-ROS point.

We elaborate on the relationship between (i) the RADA algorithmic framework and RADA-RGD in particular, which are designed for problem (1), and (ii) the Riemannian augmented Lagrangian method (RALM) [14, 83, 13] and RADMM [40], which are designed for problem (7). Roughly speaking, the subproblem for updating x in RALM is the same as that in RADA when certain regularization parameter in the aforementioned value function is zero. However, the stopping criterion for the subproblem in RADA is less stringent than that in RALM. Moreover, we show that RADA-RGD with certain parameter setting gives rise to a new Riemannian symmetric Gauss-Seidel ADMM (sGS-ADMM). These similarities and differences provide important insights into understanding the superior performance of our proposed algorithms over existing ones [83, 40].

We present numerical results on SPCA, FPCA, and SSC to demonstrate the advantages of our proposed algorithms over existing state-of-the-art algorithms for the corresponding problems. Our numerical results on SPCA show that when compared against the nested-loop algorithms ManPG [11] and RALM [83], which employ second-order methods to solve their respective subproblems, our proposed RADA-RGD can return a solution of comparable quality with much lower computational cost. In addition, our numerical results on SSC and FPCA show that our proposed algorithms outperform several state-of-the-art single-loop first-order algorithms [4, 74, 51, 40] in terms of both solution quality and speed.

1.4 Organization

The rest of the paper is organized as follows. We first introduce the notation and cover some preliminaries on Riemannian optimization and convex analysis in Section 2. We then present our proposed algorithmic framework and establish its iteration complexity in Section 3. Subsequently, we develop two customized algorithms within the proposed framework in Section 4. We reveal the connections between the proposed algorithms and some existing ones in Section 5. Then, we present numerical results on SPCA, FPCA, and SSC to illustrate the efficiency of the proposed algorithms in Section 6. Finally, we draw some conclusions in Section 7.

2 Notation and Preliminaries

We first introduce the notation and some basic concepts in Riemannian optimization [1, 9]. We use ⟨⋅,⋅⟩ and ∥⋅∥ to denote the standard inner product and its induced norm on the Euclidean space ℰ, respectively. Let 𝒳 be a subset of ℰ. We use dist⁡(y,𝒳)=infx∈𝒳‖y−x‖ to denote the distance from y∈ℰ to 𝒳 and conv​𝒳 to denote the convex hull of 𝒳. For a linear mapping 𝒜:ℰ1→ℰ2, we use 𝒜⊤:ℰ2→ℰ1 to denote its adjoint mapping. Let ℳ be a submanifold embedded in ℰ and Tx​ℳ denote the tangent space to ℳ at x∈ℳ. Throughout this paper, we take the standard inner product ⟨⋅,⋅⟩ on the Euclidean space ℰ as the Riemannian metric on ℳ. Then, for a smooth function f:ℰ→ℝ and a point x∈ℳ, the Riemannian gradient of f at x is given by grad​f​(x)=projTx​ℳ​(∇f​(x)), where proj𝒳​(⋅) is the Euclidean projection operator onto a closed set 𝒳 and ∇f​(x) is the Euclidean gradient of f at x. A retraction at x∈ℳ is a smooth mapping Rx:Tx​ℳ→ℳ satisfying (i) Rx​(𝟎x)=x, where 𝟎x is the zero element in Tx​ℳ; (ii) dd​t​Rx​(t​v)|t=0=v for all v∈Tx​ℳ.

Now, we review some basic notions in convex analysis [61, 5]. Let g:ℰ→(−∞,+∞] be a proper closed convex function whose domain is given by dom​g:={x∈ℰ∣g⁡(x)<+∞}. Its conjugate function is defined as g∗​(y):=maxx∈ℰ⁡{⟨y,x⟩−g⁡(x)}. The function g is called λ-strongly convex if g(⋅)−λ2∥⋅∥2 is convex. For a given constant λ>0, the proximal mapping and the Moreau envelope of g are defined as

proxλ​g​(x)=\argminu∈ℰ​{g⁡(u)+12​λ​‖x−u‖2}

and

Mλ​g​(x)=minu∈ℰ⁡{g⁡(u)+12​λ​‖x−u‖2},

respectively. For any x∈ℰ, we have the following so-called Moreau decomposition and Moreau envelope decomposition:

x=proxλ​g​(x)+λ​proxg∗/λ​(xλ), (8)

12​λ​‖x‖2=Mλ​g​(x)+Mg∗/λ​(xλ). (9)

The following theorem characterizes the smoothness of the Moreau envelope.

([5, Theorems 6.39 and 6.60]) Let g:ℰ→(−∞,+∞] be a proper closed convex function and λ>0. Then, for any x∈ℰ,

∇Mλ​g​(x)=1λ​(x−proxλ​g​(x))=proxg∗/λ​(xλ)∈∂g⁡(proxλ​g​(x)). (10)

Moreover, ∇Mλ​g is λ−1-Lipschitz continuous.

3 Proposed Algorithmic Framework

In this section, we propose a flexible algorithmic framework for solving the Riemannian NC-L minimax problem (1) and establish its iteration complexity.

3.1 Proposed RADA Algorithmic Framework

A widely used idea for tackling the minimax problem (1) is to alternately update the variables x and y by solving an appropriate minimization problem for x and an appropriate maximization problem for y. However, even a small perturbation in x may cause a drastic change in the solution to the problem maxy∈ℰ2⁡F⁡(x,y) [67]. Inspired in part by [53, 58, 20, 76], we address this issue by adding quadratic regularization terms to the objective function and defining the following value function at the k-th iteration:

Φk​(x):=maxy∈ℰ2⁡{F⁡(x,y)−λ2​‖y‖2−βk2​‖y−yk‖2}. (11)

Here, λ>0 is the regularization parameter, βk≥0 is the proximal parameter, and yk∈dom​h is the proximal center. Intuitively, the quadratic regularization term λ2​‖y‖2 ensures the smoothness of Φk (as shown later in (18)), while the proximal term βk2​‖y−yk‖2 stabilizes the update of y.

At the k-th iteration, our algorithmic framework computes the update xk+1 by approximately solving the subproblem

minx∈ℳ⁡Φk​(x). (12)

Specifically, we need xk+1 to satisfy the sufficient decrease condition

Φk​(xk+1)−Φk​(xk)≤−min⁡{cλ​‖grad​Φk​(xk)‖,c′}​‖grad​Φk​(xk)‖+νk, (13)

where cλ>0 is a predetermined constant on the order of λ (i.e., cλ=𝒪⁡(λ)), c′>0 is another predetermined constant, and

{νk}∈𝕊:={{νk}|∑k=1+∞νk<+∞,νk≥0,k=1, 2,…}

is a nonnegative summable sequence. Then, we compute the update

yk+1 =\argmaxy∈ℰ2​{F⁡(xk+1,y)−λ2​‖y‖2−βk2​‖y−yk‖2} (14)

=proxh/(λ+βk)​(𝒜⁡(xk+1)+βk​ykλ+βk).

The above update of y can also be regarded as performing a proximal gradient ascent step of the regularized function F(xk+1,⋅)−λ2∥⋅∥2. The proposed RADA algorithmic framework is formally presented in Algorithm 1.

Algorithm 1 RADA algorithmic framework for problem (1)

Several remarks on the proposed Algorithm 1 are in order. First, the requirement on the update βk+1 is quite flexible. A simple approach is to set βk=β1/kρ for all k≥1. Alternatively, we can adaptively update βk for improved practical efficiency. For instance, given constants τ1∈(0,1) and τ2∈(0,1), we may set

βk+1=β1(k+1)(k+1)ρwithβ1(k+1)={τ2​β1(k),if ​δk+1≥τ1​δk;β1(k),if ​δk+1<τ1​δk, (15)

where β1(1)=β1 and δk+1 is defined as

δk+1:=‖(λ+βk)​yk+1−βk​yk‖∞,∀k≥0. (16)

Here, we adopt the convention that β0=β1 and y0=y1 when computing δ1. It is simple to see that the update βk+1 in (15) satisfies 0≤βk+1≤β1/(k+1)ρ for all k≥1. As will be discussed at the end of Section 5.1, our proposed Algorithm 1 is closely related to RALM. In particular, the parameter δk+1 in (16) can be viewed as a primal residual in RALM.

Second, the sufficient decrease condition (13) is easy to satisfy, making the proposed algorithmic framework computationally efficient. For instance, we can find the desired xk+1∈ℳ by performing a simple Riemannian gradient descent step on x with a suitable stepsize. To do this, we need to establish the smoothness of Φk. By the fact that

⟨𝒜⁡(x),y⟩−λ2​‖y‖2−βk2​‖y−yk‖2+βk2​‖yk‖2

= 12​(λ+βk)​‖𝒜⁡(x)+βk​yk‖2−λ+βk2​‖y−𝒜⁡(x)+βk​ykλ+βk‖2,

we have from (11) and (9) that

Φk​(x)= f⁡(x)+maxy∈ℰ2⁡{⟨𝒜⁡(x),y⟩−h⁡(y)−λ2​‖y‖2−βk2​‖y−yk‖2} (17)

= f⁡(x)+12​(λ+βk)​‖𝒜⁡(x)+βk​yk‖2−Mh/(λ+βk)​(𝒜⁡(x)+βk​ykλ+βk)−βk2​‖yk‖2

= f⁡(x)+M(λ+βk)​h∗​(𝒜⁡(x)+βk​yk)−βk2​‖yk‖2.

Therefore, by Theorem 2.1, the value function Φk is differentiable and

∇Φk(x)=∇f(x)+∇𝒜(x)⊤proxh/(λ+βk)(𝒜⁡(x)+βk​ykλ+βk), (18)

where ∇𝒜 denotes the Jacobian mapping of 𝒜. More approaches for approximately solving the subproblem (12) are presented in Section 4. The choice of the approach mainly depends on the structure of the underlying problem. The update yk+1 in (14) is also simple since h is assumed to have a tractable proximal mapping. Therefore, our proposed algorithmic framework enjoys a low per-iteration computational cost.

Finally, we highlight the difference between our proposed RADA algorithmic framework and some related algorithms in the literature. The works [58, 53, 76, 20] for solving Euclidean NC-L minimax problems minimize certain surrogate of F⁡(⋅,y) with a fixed y to update x, i.e.,

xk+1≈\argminx∈ℳ​{F⁡(x,yk)−λ2​‖yk‖2},

where ℳ is a compact convex set. By contrast, the proposed value function (11) at x is the maximization of a regularized F⁡(x,⋅), and the subproblem (12) used to update xk+1 is itself a minimax problem, i.e.,

xk+1≈\argminx∈ℳ​{maxy∈ℰ2⁡{F⁡(x,y)−λ2​‖y‖2−βk2​‖y−yk‖2}},

which can be viewed as a surrogate of problem (1). Thus, our proposed RADA is fundamentally different from the approaches in [58, 53, 76, 20]. To provide more insights into why the former could outperform the latter, we elaborate on the difference between RADA and the ARPGDA algorithm developed in [74], which can be regarded as a Riemannian counterpart of the methods in [58, 53, 76, 20]. ARPGDA computes the update

xk+1=Rxk​(−ζk​gradx​F​(xk,yk)),

where gradx​F​(x,y) denotes the Riemannian gradient of F⁡(x,y) with respect to x. Then, it computes the update yk+1 according to (14). Similar to the derivation of ∇Φk​(x), we have ∇xF​(xk,yk)=∇Φk−1​(xk) and hence

gradx​F​(xk,yk)=grad​Φk−1​(xk).

Therefore, ARPGDA uses the Riemannian gradient of the value function Φk−1 at xk to compute the update xk+1. By contrast, our proposed framework uses the Riemannian gradient of the value function Φk at xk to compute the update xk+1. In particular, our proposed framework makes use of the most recent information about y (since Φk involves the latest iterate yk) to compute the update xk+1. This suggests that RADA could outperform ARPGDA.

3.2 Iteration Complexity

In this subsection, we study the convergence behavior of the proposed RADA algorithmic framework. Towards that end, let us first introduce the stationarity measures that we are interested in. The following definition is motivated by the stationarity measures used in Euclidean minimax problems [76, 41] and (Riemannian) nonsmooth composite problems [4, 40, 69].

Let ε>0 be some given constant.

The point (x,y)∈ℳ×dom​h is an ε-Riemannian-game-stationary (ε-RGS) point of problem (1) if there exists a constant γ>0 such that

max⁡{‖gradx​F​(x,y)‖,1γ​‖y−proxγ​h​(y+γ​𝒜​(x))‖}≤ε.

The point x∈ℳ is an ε-Riemannian-optimization-stationary (ε-ROS) point of problem (1) if there exists a point p∈ℰ2 such that

max{dist(0,gradf(x)+projTx​ℳ(∇𝒜(x)⊤∂h∗(p))),‖p−𝒜(x)‖}≤ε,

where ∂h∗ is the subdifferential of h∗.

For simplicity, denote

Fk​(x,y) :=F⁡(x,y)−λ2​‖y‖2−βk2​‖y−yk‖2, (19)

yk+12 :=\argmaxy∈ℰ2​Fk​(xk,y)=proxh/(λ+βk)​(𝒜⁡(xk)+βk​ykλ+βk), (20)

Φ⁡(x) :=maxy∈ℰ2⁡F⁡(x,y), (21)

R :=maxy∈dom​h⁡‖y‖<+∞. (22)

Recalling that ρ>1, 0≤βk≤β1/kρ, and {νk}∈𝕊 in Algorithm 1, we have

Υ:=λ​R22+∑k=1+∞(νk+2​βk​R2)<+∞. (23)

We make the following blanket assumption. {assumption} The function Φ defined in (21) is bounded from below on ℳ, namely, Φ∗:=infx∈ℳΦ⁡(x)>−∞.

Now, we are ready to establish the iteration complexity of the proposed RADA algorithmic framework for returning an ε-RGS point and an ε-ROS point of problem (1). We start with the following estimate:

Let {xk} be the sequence generated by Algorithm 1. Then, we have

Φk+1​(xk+1)−Φk​(xk)≤−min⁡{cλ​‖grad​Φk​(xk)‖,c′}​‖grad​Φk​(xk)‖+νk+2​βk​R2, (24)

where the constant cλ is given after (13).

Proof By the definition of Φk in (11) and the update formula for yk+1 in (14), we have Φk​(xk+1)=Fk​(xk+1,yk+1). Similarly, by the definition of yk+12 in (20), we have Φk+1​(xk+1)=Fk+1​(xk+1,yk+32), where we use yk+32 as a shorthand for y(k+1)+12. We now compute

Φk+1​(xk+1)−Φk​(xk+1)= Fk+1​(xk+1,yk+32)−Fk​(xk+1,yk+1)

≤(a) Fk+1​(xk+1,yk+32)−Fk​(xk+1,yk+32)

=(b) F⁡(xk+1,yk+32)−λ2​‖yk+32‖2−βk+12​‖yk+32−yk+1‖2

−F⁡(xk+1,yk+32)+λ2​‖yk+32‖2+βk2​‖yk+32−yk‖2

= βk2​‖yk+32−yk‖2−βk+12​‖yk+32−yk+1‖2​≤(c)​2​βk​R2, (25)

where (a) comes from the optimality of yk+1 in (14), (b) holds by the definition of Fk+1 in (19), and (c) is due to (22) and the fact that yk,yk+32∈dom​h. By combining (25) and (13), we obtain (24).

Based on Lemma 3.2, we can establish the iteration complexity of Algorithm 1 for returning an ε-RGS point of problem (1).

Given a constant ε>0, let {(xk,yk)} be the sequence generated by Algorithm 1 with λ=ε/(2​R). Suppose that Assumption 3.2 holds. Then, there exists a k^≤K with

K:=⌈(4​β1​Rε)1ρ⌉+⌈Φ1​(x1)−Φ∗+Υmin⁡{cλ​ε,c′}​ε⌉ (26)

such that (xk^,yk^+12) is an ε-RGS point of problem (1).

Proof Given any γ>0, we first show that yk+12 satisfies

1γ​‖yk+12−proxγ​h​(yk+12+γ​𝒜​(xk))‖≤ε,∀k≥K1:=⌈(4​β1​Rε)1ρ⌉. (27)

For brevity, denote yk+12′=proxγ​h​(yk+12+γ​𝒜​(xk)). It is easy to see that

yk+12′=\argmaxy∈ℰ2​{F⁡(xk,y)−12​γ​‖y−yk+12‖2}. (28)

Due to the optimality of yk+12 in (20), we have

0≤Fk​(xk,yk+12)−Fk​(xk,yk+12′). (29)

Since F(xk,⋅)−12​γ∥⋅−yk+12∥2 is γ−1-strongly concave, it follows from the optimality of yk+12′ in (28) that

12​γ​‖yk+12−yk+12′‖2≤F⁡(xk,yk+12′)−12​γ​‖yk+12′−yk+12‖2−F⁡(xk,yk+12).

This, together with (29), further implies that

1γ​‖yk+12−yk+12′‖2≤ (Fk​(xk,yk+12)−F⁡(xk,yk+12))−(Fk​(xk,yk+12′)−F⁡(xk,yk+12′))

=(a) λ2​(‖yk+12′‖2−‖yk+12‖2)+βk2​(‖yk+12′−yk‖2−‖yk+12−yk‖2)

= λ2​⟨yk+12′+yk+12,yk+12′−yk+12⟩+βk2​⟨yk+12′+yk+12−2​yk,yk+12′−yk+12⟩

≤ (λ2​‖yk+12′+yk+12‖+βk2​‖yk+12′+yk+12−2​yk‖)​‖yk+12−yk+12′‖,

where (a) uses the definition of Fk in (19). Cancelling the term ‖yk+12−yk+12′‖ on the two ends of the above inequalities yields

1γ​‖yk+12−yk+12′‖≤λ2​‖yk+12′+yk+12‖+βk2|yk+12′+yk+12−2​yk|≤(a)​λ​R+2​βk​R​≤(b)​ε2+ε2=ε,

where (a) is due to (22) and the fact that yk,yk+12,yk+12′∈dom​h, and (b) uses λ=ε/(2​R), 0≤βk≤β1/kρ, and k≥K1=⌈(4​β1​R/ε)1ρ⌉.

Next, we prove the existence of some k^∈{K1,K1+1,…,K} such that xk^ satisfies

‖gradx​F​(xk^,yk^+12)‖≤ε. (30)

Note that for any x∈ℳ, we deduce from (11) and the expression of F⁡(x,y) in (1) that

Φk​(x) =maxy∈ℰ2⁡{f⁡(x)+⟨𝒜⁡(x),y⟩−h⁡(y)−λ2​‖y‖2−βk2​‖y−yk‖2} (31)

≥maxy∈ℰ2⁡{f⁡(x)+⟨𝒜⁡(x),y⟩−h⁡(y)−λ2​R2−2​βk​R2}

=Φ⁡(x)−λ2​R2−2​βk​R2,

where the inequality is due to (22) and the last equality follows from the definition of Φ in (21). By (24), we know that

min⁡{cλ​‖grad​Φk​(xk)‖,c′}​‖grad​Φk​(xk)‖≤Φk​(xk)−Φk+1​(xk+1)+νk+2​βk​R2.

Summing the above inequality over k=K1,K1+1,…,K yields

∑k=K1Kmin⁡{cλ​‖grad​Φk​(xk)‖,c′}​‖grad​Φk​(xk)‖

≤ Φ1​(x1)−ΦK+1​(xK+1)+∑k=1K(νk+2​βk​R2)≤Φ1​(x1)−Φ∗+Υ,

where the last inequality is due to Assumption 3.2, (23), and (31). Hence, there exists some k^∈{K1,K1+1,…,K} such that

min⁡{cλ​‖grad​Φk^​(xk^)‖,c′}​‖grad​Φk^​(xk^)‖≤Φ1​(x1)−Φ∗+ΥK−K1. (32)

Substituting (26) into (32), we have

‖grad​Φk^​(xk^)‖≤ε. (33)

The definition of yk+12 in (20), together with (18), implies that ∇Φk​(xk)=∇xF​(xk,yk+12) and thus

grad​Φk​(xk)=gradx​F​(xk,yk+12). (34)

Therefore, we see from (33) and (34) that (30) holds. Together with (27), we conclude that (xk^,yk^+12) is an ε-RGS point of problem (1).

According to Theorem 3.4, since λ=ε/(2​R) and cλ=𝒪⁡(λ) by our choice, we see that the iteration complexity of Algorithm 1 for obtaining an ε-RGS point of problem (1) is 𝒪⁡(ε−3). It is worth mentioning that our iteration complexity result in Theorem 3.4 for the RADA algorithmic framework matches not only that for the less general ARPGDA in [74], but also matches the best-known complexity results for algorithms that tackle NC-L minimax problems with x being constrained on a nonempty compact convex set instead of a Riemannian manifold [58, 65, 20]. Furthermore, our iteration complexity analysis is based on the sufficient decrease condition (13) of the proposed value function Φk. By contrast, the iteration complexity analyses in [58, 76, 74] rely on the descent properties of the function F⁡(⋅,yk) and a different potential function, which do not extend to RADA.

Next, we establish the iteration complexity of Algorithm 1 for returning an ε-ROS point of problem (1).

Consider the setting of Theorem 3.4. Then, xk^ is an ε-ROS point of problem (1), where k^ is the index given in Theorem 3.4.

Proof The proof is similar to that of Theorem 3.4. From (33), we know there exists some k^∈{K1,K1+1,…,K} such that

‖grad​Φk^​(xk^)‖≤ε. (35)

Next, define pk=prox(λ+βk)​h∗​(𝒜⁡(xk)+βk​yk). From (10), (18), and (20), we have

grad​Φk​(xk) =gradf(xk)+projTxk​ℳ(∇𝒜(xk)⊤yk+12)andyk+12∈∂h∗(pk). (36)

Then, it follows from (35) and (36) that

dist(0,gradf(xk^)+projTxk^​ℳ(∇𝒜(xk^)⊤∂h∗(pk^))) (37)

≤ ‖gradf(xk^)+projTx​ℳ(∇𝒜(xk^)⊤yk^+12)‖=‖gradΦk^(xk^)‖≤ε.

On the other hand, by the Moreau decomposition (8), we have

‖pk^−𝒜⁡(xk^)‖= ‖prox(λ+βk^)​h∗​(𝒜⁡(xk^)+βk^​yk^)−𝒜⁡(xk^)‖ (38)

= ‖βk^​yk^−(λ+βk^)​proxh/(λ+βk^)​(𝒜⁡(xk^)+βk^​yk^λ+βk^)‖

≤(a) (λ+2​βk^)​R​≤(b)​ε,

where (a) uses proxh/(λ+βk^)​(⋅)∈dom​h and (22), and (b) uses λ=ε/(2​R), 0≤βk^≤β1/k^ρ, and k^≥K1=⌈(4​β1​R/ε)1ρ⌉. Combining (37) and (38), we conclude that xk^ is an ε-ROS point of problem (1).

Theorem 3.6 shows that the iteration complexity of Algorithm 1 for obtaining an ε-ROS point of problem (1) is 𝒪⁡(ε−3) since λ=ε/(2​R) and cλ=𝒪⁡(λ). This matches the best-known complexity results for the DSGM in [4] and the ManIAL in [13]. However, the results in [4, 13] apply only to the case where 𝒜 is a linear mapping, whereas our results work for a general nonlinear mapping 𝒜.

Before we leave this section, let us point out that by Theorems 3.4 and 3.6, Algorithm 1 can find both an ε-RGS point and an ε-ROS point of problem (1) within 𝒪⁡(ε−3) iterations. This is noteworthy, as the results in [46, 77] for the Euclidean setting suggest that in general extra computation is needed to convert between an ε-RGS point and an ε-ROS point. Moreover, existing works on NC-L minimax problems [58, 20, 74] only establish the iteration complexity for obtaining an ε-game-stationary point. The iteration complexity for obtaining an ε-optimization-stationary point is not addressed in these works.

4 Two Efficient Algorithms for Problem (1)

In this section, we propose two efficient algorithms within the proposed RADA algorithmic framework, called RADA-PGD and RADA-RGD in Algorithms 2 and 3, respectively, for solving the Riemannian NC-L minimax problem (1). These two algorithms differ in the first-order method used for approximately solving the subproblem (12) to find a point xk+1∈ℳ that satisfies the sufficient decrease condition (13).

Algorithm 2 RADA-PGD for problem (1)

When the projection operator onto a closed manifold ℳ is tractable, which is the case in the motivating applications introduced in Section 1.1, we can employ the nonconvex projected gradient method to compute xk+1 that satisfies the sufficient decrease condition (13). Noting that Φk is differentiable as shown in (18), it is straightforward to perform an ordinary projected gradient step with a suitable constant stepsize depending on the Lipschitz constant of ∇Φk. It is also possible to perform multiple projected gradient steps to achieve a more accurate solution to (12) while satisfying (13). Within the RADA algorithmic framework, the first customized algorithm is RADA with projected gradient descent (RADA-PGD), which is formally presented in Algorithm 2. Here, ℤ+ is the set of all positive integers, and the integer Tk∈[1,T¯] denotes the number of the projected gradient descent steps to be performed at the k-th iteration with T¯ being a predetermined positive integer. The stepsize ζk,t is set as ℓk−1, where ℓk is to be specified later in (46).

An alternative approach is to perform a Riemannian gradient descent step to obtain xk+1 satisfying the sufficient decrease condition (13). For this method, determining a suitable constant stepsize is not as easy as in RADA-PGD. Hence, for practical efficiency, we utilize the Barzilai-Borwein (BB) stepsize [3, 71, 30, 24, 29, 31, 16] as the initial stepsize and perform a backtracking line search to find a stepsize satisfying (13). Similar to Algorithm 2, we can also perform multiple Riemannian gradient steps to achieve a more accurate solution to (12). The corresponding algorithm is named RADA with Riemannian gradient descent (RADA-RGD) and is presented in Algorithm 3, where the integer Tk∈[1,T¯] denotes the number of the Riemannian gradient descent steps to be performed at the k-th iteration and T¯ is a predetermined positive integer. The Riemannian BB stepsize ζk,t+1BB in Algorithm 3 is set as

ζk,t+1BB={‖xk,t+1−xk,t‖2|⟨xk,t+1−xk,t,vk,t⟩|, for odd ​t;|⟨xk,t+1−xk,t,vk,t⟩|‖vk,t‖2, for even ​t, (39)

where vk,t=grad​Φk​(xk,t+1)−grad​Φk​(xk,t).

Note that both RADA-PGD and RADA-RGD are single-loop algorithms, as the inner iteration number Tk is bounded by a predetermined positive integer T¯. The choice of Tk is also flexible and can be adapted to the specific problem structure. Empirically, if the value function Φk is benign, for example, if the Lipschitz constant of its gradient can be well estimated, then setting Tk≡1 is often sufficient. Otherwise, we suggest setting a slightly larger Tk to ensure a more substantial decrease in the value of Φk, thereby improving the practical performance of RADA-PGD and RADA-RGD. This is particularly important because the subproblem (12) is itself a minimax problem and can be viewed as a surrogate of problem (1). Thanks to the proposed value function Φk and the flexible setting of Tk, our approaches to updating xk+1 are expected to be more efficient than simply performing a gradient descent step on certain surrogate of F⁡(⋅,yk), as is done in [76, 58, 53, 20, 74]. This will be demonstrated in Section 6. Moreover, even when adapted to the Euclidean setting, our proposed RADA algorithmic framework and the associated RADA-PGD and RADA-RGD are new.

Algorithm 3 RADA-RGD for problem (1)

Φk​(Rxk,t​(−ζk,t​ηj​grad​Φk​(xk,t)))−Φk​(xk,t)≤−c1​ζk,t​ηj​‖grad​Φk​(xk,t)‖2+νkTk, (40)

ζk,t+1=min⁡{max⁡{ζk,t+1BB,ζmin},ζmax/‖grad​Φk​(xk,t+1)‖}. (41)

4.1 Verification of the Sufficient Decrease Condition (13)

In this subsection, we show that the update xk+1 in RADA-PGD and RADA-RGD satisfies the sufficient decrease condition (13). Therefore, the iteration complexity of these two algorithms can be directly obtained from the results in Section 3.2.

First, we make the following (somewhat standard) assumptions for our analysis; see, e.g., [12, 11, 8, 27, 28, 48, 15]. {assumption} The level set

Ωx1:={x∈ℳ∣Φ⁡(x)≤Φ⁡(x1)+Υ}

is compact, where Φ and Υ are defined in (21) and (23), respectively, and x1 is the initial point of Algorithm 1. {assumption} The function f and the mapping 𝒜 in (1) satisfy the following conditions:

The function f:ℰ1→ℝ is continuously differentiable and satisfies the descent property over ℳ, i.e.,

f⁡(x′)≤f⁡(x)+⟨∇f​(x),x′−x⟩+Lf2​‖x′−x‖2,∀x,x′∈ℳ. (42)

The mapping 𝒜:ℰ1→ℰ2 and its Jacobian mapping ∇𝒜 are L𝒜0-Lipschitz and L𝒜1-Lipschitz over conv​ℳ, respectively. In other words, for any x,x′∈conv​ℳ, we have

‖𝒜⁡(x)−𝒜⁡(x′)‖≤L𝒜0​‖x−x′‖, (43a)

‖∇𝒜​(x)−∇𝒜​(x′)‖≤L𝒜1​‖x−x′‖. (43b)

Moreover, the mapping ∇𝒜 is bounded over conv​ℳ, i.e.,

ρ𝒜:=supx∈conv​ℳ‖∇𝒜​(x)‖<+∞. (44)

In our subsequent development, we allow for the possibility that the retraction at some x∈ℳ is defined only locally; see [8, Definition 2.1]. As such, we make the following assumption, which stipulates that the retraction is well defined on a compact subset of the tangent bundle; cf. [8, Remark 2.2]. {assumption} There exists a constant ϱ>0 such that for any (x,v)∈𝒱:={(x,v)∣x∈Ωx1,v∈Txℳ,∥v∥≤ϱ}, the point Rx​(v)∈ℳ is well defined. The following lemma, which is extracted from [8, Appendix B], shows that the retraction satisfies the first- and second-order boundedness conditions on 𝒱.

Suppose that Assumptions 4.1 and 4.1 hold. Then, there exist constants α1,α2>0 such that

‖Rx​(v)−x‖≤α1​‖v‖and‖Rx​(v)−x−v‖≤α2​‖v‖2

for any (x,v)∈𝒱.

Based on Lemma 4.1, we can establish the descent property of Φk as follows.

Suppose that Assumptions 4.1, 4.1, and 4.1 hold. Then, the value function Φk in (11) satisfies the following properties:

Euclidean Descent. For any x,x′∈ℳ, we have

Φk​(x′)≤Φk​(x)+⟨∇Φk​(x),x′−x⟩+ℓk2​‖x′−x‖2, (45)

where

ℓk=Lf+R​L𝒜1+ρ𝒜​L𝒜0λ+βk. (46)

Riemannian Descent. For any (x,v)∈𝒱, we have

Φk​(Rx​(v))≤Φk​(x)+⟨grad​Φk​(x),v⟩+Lk​(x)2​‖v‖2, (47)

where

Lk​(x)=ℓk​α12+2​(‖∇f​(x)‖+ρ𝒜​R)​α2. (48)

Proof We first prove (45). Let Pk​(x)=Φk​(x)−f⁡(x). By (18), we have

∇Pk(x)=∇𝒜(x)⊤ℬ(x),

where

ℬ⁡(x)=proxh/(λ+βk)​(𝒜⁡(x)+βk​ykλ+βk).

We claim that ∇Pk is Lipschitz continuous. Indeed, for any x,x′∈conv​ℳ, we have

‖∇Pk​(x)−∇Pk​(x′)‖= ∥∇𝒜(x)⊤ℬ(x)−∇𝒜(x′)⊤ℬ(x′)∥

≤ ∥(∇𝒜(x)−∇𝒜(x′))⊤ℬ(x)∥+∥∇𝒜(x′)⊤(ℬ(x)−ℬ(x′))∥

≤(a) L𝒜1​‖ℬ⁡(x)‖⋅‖x−x′‖+ρ𝒜​‖ℬ⁡(x)−ℬ⁡(x′)‖

≤(b) L𝒜1​‖ℬ⁡(x)‖⋅‖x−x′‖+ρ𝒜λ+βk​‖𝒜⁡(x)−𝒜⁡(x′)‖

≤(c) (R​L𝒜1+ρ𝒜​L𝒜0λ+βk)​‖x−x′‖,

where (a) is due to (43b) and (44), (b) is due to the nonexpansiveness of the proximal operator, and (c) is due to ℬ⁡(x)∈dom​h, (22), and (43a). Therefore, by [55, Lemma 1.2.3], we know that Pk satisfies the Euclidean descent property, namely, for any x,x′∈ℳ,

Pk​(x′)≤Pk​(x)+⟨∇Pk​(x),x′−x⟩+12​(R​L𝒜1+ρ𝒜​L𝒜0λ+βk)​‖x′−x‖2. (49)

Combining (49) and (42) and noting that Φk=f+Pk, we conclude that Φk satisfies the Euclidean descent property (45).

Moreover, using Lemma 4.1 and following a similar analysis as that in [8, Appendix B], we know that Φk satisfies the Riemannian descent property (47). This completes the proof.

To prove that xk+1 in Algorithm 2 satisfies the sufficient decrease condition (13), we need the following lemma.

Let {xk,t}t=1Tk+1 be the sequence generated by Algorithm 2 with ζk,t=ℓk−1 at the k-th iteration. Suppose that Assumptions 4.1, 4.1, and 4.1 hold. If xk,t∈Ωx1, then

Φk​(xk,t+1)−Φk​(xk,t)≤−min⁡{‖grad​Φk​(xk,t)‖2​Lk​(xk,t),ϱ2}​‖grad​Φk​(xk,t)‖.

Proof By the definition of xk,t+1 in Algorithm 2, we have

xk,t+1∈ \argminx∈ℳ‖x−xk,t+1ℓk∇Φk(xk,t)‖2 (50)

= \argminx∈ℳ​{⟨∇Φk​(xk,t),x−xk,t⟩+ℓk2​‖x−xk,t‖2}.

Consider xk,t+1′=Rxk,t​(dk,t) with dk,t=−ck,t​grad​Φk​(xk,t), where

ck,t:=min⁡{1Lk​(xk,t),ϱ‖grad​Φk​(xk,t)‖}. (51)

We know that (xk,t,dk,t)∈𝒱 from (51) where 𝒱 is defined in Assumption 4.1. By the optimality of xk,t+1 in (50), we get

⟨∇Φk​(xk,t),xk,t+1−xk,t⟩+ℓk2​‖xk,t+1−xk,t‖2

≤ ⟨∇Φk​(xk,t),xk,t+1′−xk,t⟩+ℓk2​‖xk,t+1′−xk,t‖2

= ⟨∇Φk​(xk,t),dk,t⟩+⟨∇Φk​(xk,t),xk,t+1′−xk,t−dk,t⟩+ℓk2​‖xk,t+1′−xk,t‖2

≤ ⟨∇Φk​(xk,t),dk,t⟩+‖∇Φk​(xk,t)‖⋅‖xk,t+1′−xk,t−dk,t‖+ℓk2​‖xk,t+1′−xk,t‖2

≤(a) −ck,t​‖grad​Φk​(xk,t)‖2+α2​(‖∇f​(x)‖+ρ𝒜​R)​‖dk,t‖2+α12​ℓk2​‖dk,t‖2

=(b) −ck,t​(1−Lk​(xk,t)2​ck,t)​‖grad​Φk​(xk,t)‖2,

where (a) follows from Lemma 4.1, proxh/(λ+βk)​((𝒜⁡(x)+βk​yk)/(λ+βk))∈dom​h, (18), (22), and (44), and (b) is due to (48). In addition, applying (45) with x′=xk,t+1 and x=xk,t yields

Φk​(xk,t+1)−Φk​(xk,t)≤⟨∇Φk​(xk,t),xk,t+1−xk,t⟩+ℓk2​‖xk,t+1−xk,t‖2

≤ −ck,t​(1−Lk​(xk,t)2​ck,t)​‖grad​Φk​(xk,t)‖2​≤(a)−ck,t2​‖grad​Φk​(xk,t)‖2,

where (a) is due to (51). The desired result follows directly from substituting (51) into the above inequalities.

Based on Lemma 4.4, we now present Proposition 4.6, which shows that Algorithm 2 is within the RADA algorithmic framework.

Let {xk} be the sequence generated by Algorithm 2 with ζk,t=ℓk−1. Suppose that Assumptions 4.1, 4.1, and 4.1 hold. Then, the sequence {xk} lies in the compact level set Ωx1. Moreover, for k=1,2,…, we have

Φk​(xk+1)−Φk​(xk)≤−min⁡{‖grad​Φk​(xk)‖2​Lk​(xk),ϱ2}​‖grad​Φk​(xk)‖. (52)

Proof We apply mathematical induction. First, we show that (52) holds when k=1. Since x1,1=x1∈Ωx1, using Lemma 4.4 with t=1 and the definition of ck,t in (51), we have

Φ1​(x1,2)−Φ1​(x1,1)≤−min⁡{‖grad​Φ1​(x1,1)‖2​L1​(x1,1),ϱ2}​‖grad​Φ1​(x1,1)‖≤0.

This gives

Φ⁡(x1,2)​≤(a)​Φ1​(x1,2)+(λ2+2​β1)​R2​≤(b)​Φ1​(x1,1)+Υ​≤(c)​Φ​(x1)+Υ, (53)

where (a) follows from (31), (b) uses the definition of Υ in (23), and (c) holds by the definition of Φ in (21) and the fact that

Φ1​(x)= maxy∈ℰ2⁡{f⁡(x)+⟨𝒜⁡(x),y⟩−h⁡(y)−λ2​‖y‖2−β12​‖y−y1‖2} (54)

≤ maxy∈ℰ2⁡{f⁡(x)+⟨𝒜⁡(x),y⟩−h⁡(y)}=Φ⁡(x).

Hence, we have x1,2∈Ωx1 from (53). Repeating the above procedure, we further know that

Φ1​(x1,t+1)−Φ1​(x1,t)≤−min⁡{‖grad​Φ1​(x1,t)‖2​L1​(x1,t),ϱ2}​‖grad​Φ1​(x1,t)‖≤0 (55)

and x1,t+1∈Ωx1 for t=1,2,…,T1. Summing the above inequality over t=1,2,…,T1 yields

Φ1​(x1,T1+1)−Φ1​(x1,1)= ∑t=1T1(Φ1​(x1,t+1)−Φ1​(x1,t)) (56)

≤ −min⁡{‖grad​Φ1​(x1)‖2​L1​(x1),ϱ2}​‖grad​Φ1​(x1)‖.

Recalling that x1,1=x1 and x1,T1+1=x2, it follows from (56) that (52) holds for k=1.

Next, we assume that xj∈Ωx1 for some j≥2 and (52) holds for k=1,2,…,j−1. It remains to prove that xj+1∈Ωx1 and (52) holds for k=j. Combining (25) and (52) with k=1,2,…,j−1 yields

Φk+1(xk+1)−Φk(xk)≤2βkR2,∀k=1,2,…,j−1. (57)

Since xj,1=xj∈Ωx1, Lemma 4.4 implies that

Φj​(xj,2)≤Φj​(xj,1)=Φj​(xj). (58)

Similar to the proof for j=1, we have

Φ⁡(xj,2)≤ Φj​(xj,2)+(λ2+2​βj)​R2 (59)

≤(a) Φj​(xj)+(λ2+2​βj)​R2

= Φ1​(x1)+∑k=2j(Φk​(xk)−Φk−1​(xk−1))+(λ2+2​βj)​R2

≤(b) Φ1​(x1)+λ2​R2+∑k=1+∞(νk+2​βk​R2)​≤(c)​Φ​(x1)+Υ,

where (a) is from (58), (b) is due to (57) and βj≤β1, and (c) is due to (54) and the definition of Υ in (23). Hence, we get xj,2∈Ωx1. Similar to (55) and (56), we deduce from Lemma 4.4 that

Φj​(xj,t+1)−Φj​(xj,t)≤−min⁡{‖grad​Φj​(xj,t)‖2​Lj​(xj,t),ϱ2}​‖grad​Φj​(xj,t)‖≤0 (60)

for t=1, 2,…,Tj. It follows from xj,1=xj and xj,Tj+1=xj+1 that

Φj​(xj+1)−Φj​(xj)=∑t=1Tj(Φj​(xj,t+1)−Φj​(xj,t))≤−min⁡{‖grad​Φj​(xj)‖2​Lj​(xj),ϱ2}​‖grad​Φj​(xj)‖.

Moreover, similar to (59), we have

Φ⁡(xj,Tj+1)≤ Φj​(xj,Tj+1)+(λ2+2​βj)​R2

≤(a) Φj​(xj)+(λ2+2​βj)​R2

= Φ1​(x1)+∑k=2j(Φk​(xk)−Φk−1​(xk−1))+(λ2+2​βj)​R2

≤(b) Φ1​(x1)+λ2​R2+∑k=1+∞(νk+2​βk​R2)​≤(c)​Φ​(x1)+Υ,

where (a) is from (60) and xj,1=xj, (b) is due to (57) and βj≤β1, and (c) is due to (54) and the definition of Υ in (23). Therefore, noting that xj+1=xj,Tj+1, we have xj+1∈Ωx1. Moreover, similar to the derivation of (56), we know from (60) that (52) holds for k=j. This completes the inductive step and hence also the proof.

Denote

ρΦ:=maxx∈Ωx1⁡‖∇f​(x)‖+ρ𝒜​RandL¯:=ℓ1​α12+2​ρ𝒜​R​α2, (61)

where ρ𝒜 and R are defined in (44) and (22), respectively. Then, for any x∈Ωx1 and k=1, 2,…, we have

‖grad​Φk​(x)‖≤‖∇Φk​(x)‖≤ρΦ (62)

from (18), (22), and (44). Moreover, we have

0<L¯≤Lk​(x) (63)

by substituting 0≤βk≤β1 into (48). To prove that xk+1 in Algorithm 3 satisfies the sufficient decrease condition (13), we need the following lemma.

Let {xk,t}t=1Tk+1 be the sequence generated by Algorithm 3 with ζmax≤ϱ at the k-th iteration, where ϱ is given in Assumption 4.1. Suppose that Assumptions 4.1, 4.1, and 4.1 hold. If xk,t∈Ωx1, then

Φk​(xk,t+1)−Φk​(xk,t)≤−c1​min⁡{ρΦ−1​ζmax​L¯,ζmin​L¯,2​η​(1−c1)}Lk​(xk,t)​‖grad​Φk​(xk,t)‖2+νkTk. (64)

Proof Let j^ be the smallest nonnegative integer satisfying (40) and ζ^k,t:=ζk,t​ηj^. First, if j^≥1, then by the backtracking line search procedure, we have

Φk​(xk,t+)−Φk​(xk,t)>−c1​ζ^k,tη​‖grad​Φk​(xk,t)‖2+νkTk≥−c1​ζ^k,tη​‖grad​Φk​(xk,t)‖2, (65)

where xk,t+:=Rxk,t​(−ζ^k,t​η−1​grad​Φk​(xk,t)). Since ζ^k,t​η−1=ζk,t​ηj^−1≤ζk,t and (41) implies that

ζk,t​‖grad​Φk​(xk,t)‖≤ζmax≤ϱ,

we get (xk,t,−ζ^k,t​η−1​grad​Φk​(xk,t))∈𝒱, where 𝒱 is defined in Assumption 4.1. By Lemma 4.2 (ii), we have

Φk​(xk,t+)−Φk​(xk,t)≤−(1−ζ^k,t​Lk​(xk,t)2​η)​ζ^k,tη​‖grad​Φk​(xk,t)‖2. (66)

Combining (65) and (66) yields

ζ^k,t≥2​η​(1−c1)Lk​(xk,t).

Second, if j^=0, then ζ^k,t=ζk,t. We bound

ζ^k,t≥min⁡{ζk,t,2​η​(1−c1)Lk​(xk,t)}≥min⁡{ρΦ−1​ζmax​L¯,ζmin​L¯,2​η​(1−c1)}Lk​(xk,t), (67)

where the second inequality holds by (63), (62), and ζk,t≥min⁡{ζmin,ζmax/‖grad​Φk​(xk,t+1)‖} as defined in (41). By substituting (67) into (40), we get (64), as desired.

Based on Lemma 4.8, we have the following proposition, which shows that Algorithm 3 is within the RADA algorithmic framework.

Let {xk} be the sequence generated by Algorithm 3. Suppose that Assumptions 4.1, 4.1, and 4.1 hold. Then, the sequence {xk} lies in the compact level set Ωx1. Moreover, for k=1,2,…, we have

Φk​(xk+1)−Φk​(xk)≤−c1​min⁡{ρΦ−1​ζmax​L¯,ζmin​L¯,2​η​(1−c1)}Lk​(xk,t),‖grad​Φk​(xk)‖2+νk.

Proof The result is obtained by using the same techniques as in the proof of Proposition 4.6.

Recalling that βk≥0, we know from (46), (48), and (61) that

Lk​(xk)≤L¯:=α12​(Lf+R​L𝒜1+ρ𝒜​L𝒜0λ)+2​ρΦ​α2=𝒪⁡(λ−1).

Therefore, from Propositions 4.6 and 4.10, we see that the sufficient decrease condition (13) holds with

cλ=12​L¯andc′=ϱ2

for Algorithm 2 and with

cλ=min⁡{ρΦ−1​ζmax​L¯,ζmin​L¯,2​η​(1−c1)}L¯,andc′=cλ​ρΦ

for Algorithm 3. Here, we use ‖grad​Φk​(x)‖=min⁡{‖grad​Φk​(x)‖,ρΦ} as shown by (62). Consequently, by invoking Theorems 3.4 and 3.6, we obtain the following complexity result:

Given a constant ε>0, let {(xk,yk)} be the sequence generated by Algorithms 2 or 3 with λ=ε/(2​R). Suppose that Assumptions 3.2, 4.1, 4.1, and 4.1 hold. Then, there exists an index k^≤K with K=𝒪⁡(ε−3) such that (xk^,yk^+12) is an ε-RGS point and xk^ is an ε-ROS point of problem (1).

5 Connections with Existing Algorithms

Recall that the minimax problem (1) is equivalent to the nonsmooth composite problem (7). As it turns out, the proposed RADA algorithmic framework and the proposed RADA-RGD algorithm for solving problem (1) have interesting connections to the RALMs [14, 83, 13] and RADMM [40] for solving problem (7), respectively. In what follows, we elaborate on these connections and explain why RADA and RADA-RGD can be more advantageous than RALM and RADMM for tackling problem (1).

5.1 Connection between RADA and RALM

In this subsection, we show how the updates of the proposed RADA algorithmic framework for solving the minimax problem (1) relates to those of the RALMs [14, 83, 13] for solving the equivalent nonsmooth composite problem (7).

On one hand, using (17), the subproblem (12) for obtaining the update xk+1 in RADA can be written as

minx∈ℳ⁡{f⁡(x)+M(λ+βk)​h∗​(𝒜⁡(x)+βk​yk)}. (68)

Moreover, by the Moreau decomposition (8), the update formula (14) for yk+1 in RADA can be written as

yk+1 =1λ+βk​(βk​yk+𝒜⁡(xk+1)−prox(λ+βk)​h∗​(𝒜⁡(xk+1)+βk​yk)). (69)

On the other hand, consider applying the RALM in [14, 83, 13] to the following equivalent reformulation of (7):

minx∈ℳ,q∈ℰ2⁡{f⁡(x)+h∗​(q)}s.t.𝒜⁡(x)−q=0. (70)

Here, q is an auxiliary variable. The augmented Lagrangian function associated with problem (70) is given by

ℒσ​(x,q,y):=f⁡(x)+h∗​(q)+⟨y,𝒜⁡(x)−q⟩+σ2​‖𝒜⁡(x)−q‖2,

where y is the Lagrange multiplier corresponding to the constraint 𝒜⁡(x)−q=0 and σ>0 is the penalty parameter. At the k-th iteration, instead of utilizing the subproblem

minx∈ℳ,q∈ℰ2⁡ℒσk​(x,q,yk), (71)

RALM utilizes the equivalent subproblem

minx∈ℳ{ℒk(x):=f(x)+Mh∗/σk(𝒜(x)+ykσk)} (72)

to obtain the update xk+1, where the equivalence follows by eliminating the variable q in (71). Specifically, RALM first approximately solves (72) to obtain a point xk+1∈ℳ that satisfies

‖grad​ℒk​(xk+1)‖≤εk, (73)

where the sequence of tolerances {εk} converges to 0. Then, it updates

yk+1=yk+σk​(𝒜⁡(xk+1)−proxh∗/σk​(𝒜⁡(xk+1)+ykσk)). (74)

Now, let us compare the subproblems (68) and (72) for updating x and the formulas (69) and (74) for updating y. When λ=0 and βk=σk−1, the subproblem (68) in RADA coincides with the subproblem (72) in RALM, and the update formula (69) in RADA reduces to the update formula (74) in RALM. However, when λ>0, both (68) and (69) in RADA differ from their respective counterparts (72) and (74) in RALM. More importantly, RADA and RALM differ in their stopping criteria for the subproblems. Indeed, RADA only needs to find a solution xk+1 that satisfies the sufficient decrease condition (13), whereas RALM needs to find a solution xk+1 that satisfies the inexactness condition (73). The former offers significantly more flexibility when solving the subproblem. It should be mentioned that (68) is also related to the so-called dampened or perturbed augmented Lagrangian function (see, e.g., in [17, 37]), which is used in the design of primal-dual algorithms for linearly constrained composite optimization problems in the Euclidean space. Furthermore, when βk≡0 and λ>0, the subproblem (68) becomes

minx∈ℳ⁡{f⁡(x)+Mλ​h∗​(𝒜⁡(x))}, (75)

which corresponds to the subproblem in the Riemannian smoothing methods in [59, 4, 80].

The above relationship between the proposed RADA algorithmic framework and RALM provides important insights into the update formula (15) for βk+1 based on the parameter δk+1 defined in (16), as well as the numerical performance of the two approaches. Specifically, at the k-th iteration of RALM, we can express the update qk+1 as

qk+1=proxh∗/σk​(𝒜⁡(xk+1)+ykσk).

Using (69), we see that our Algorithm 1 has an analogous update, namely,

pk+1=prox(λ+βk)​h∗​(𝒜⁡(xk+1)+βk​yk).

Now, using (16), we compute

δk+1= ‖(λ+βk)​yk+1−βk​yk‖∞

=(a) ‖−βk​yk+(λ+βk)​proxh/(λ+βk)​(𝒜⁡(xk+1)+βk​ykλ+βk)‖∞

=(b) ‖𝒜⁡(xk+1)−prox(λ+βk)​h∗​(𝒜⁡(xk+1)+βk​yk)‖∞

= ‖𝒜⁡(xk+1)−pk+1‖∞,

where (a) is due to (14) and (b) follows from the Moreau decomposition (8). This shows that the parameter δk+1 in (16) can be viewed as a primal residual in RALM.

5.2 Connection between RADA-RGD and RADMM

In this subsection, we first elaborate on the relationship between RADA-RGD with Tk≡1 and the RADMM in [40]. We then show that RADA-RGD with Tk≡1 is equivalent to a new Riemannian sGS-ADMM.

Consider applying the RADMM in [40] to the following equivalent formulation of the smoothed problem (75):

minx∈ℳ⁡{f⁡(x)+Mλ​h∗​(p)}s.t.𝒜⁡(x)−p=0. (76)

The augmented Lagrangian function associated with problem (76) is given by

ℒ~σ​(x,p,y):=f⁡(x)+Mλ​h∗​(p)+⟨y,𝒜⁡(x)−p⟩+σ2​‖𝒜⁡(x)−p‖2,

where y is the Lagrange multiplier corresponding to the constraint 𝒜⁡(x)−p=0 and σ>0 is the penalty parameter. At the k-th iteration, RADMM computes the updates

xk+1= Rxk​(−ζk​gradx​ℒ~σk​(xk,pk,yk)), (77a)

pk+1= \argminp∈ℰ2​ℒ~σk​(xk+1,p,yk), (77b)

yk+1= yk+σk​(𝒜⁡(xk+1)−pk+1), (77c)

where the stepsize ζk>0 in (77a) and the penalty parameter σk>0 are both taken to be constants in [40].

Next, we modify RADA-RGD in Algorithm 3 by introducing an additional variable p, so that the resulting algorithm, which is shown in Algorithm 4, can be used to solve problem (7). Here, for simplicity of presentation, we omit the line search procedure in RADA-RGD.

Algorithm 4 RADA-RGD reformulation for problem (7)

pk,t =\argminp∈ℰ2​ℒ~1/βk​(xk,t,p,yk),

xk,t+1 =Rxk,t​(−ζk,t​gradx​ℒ~1/βk​(xk,t,pk,t,yk)). (78)

yk+1=yk+1βk​(𝒜⁡(xk+1)−pk+1). (79)

The key difference between (77) and Algorithm 4 with Tk≡1 is that the latter updates x and p in a symmetric Gauss-Seidel fashion, i.e., p→x→p. Therefore, Algorithm 4 with Tk≡1 can be viewed as a Riemannian sGS-ADMM (cf. [10, 43, 44] for the sGS-ADMM in Euclidean convex settings) for solving the Riemannian nonsmooth composite problem (7). To the best of our knowledge, this is the first time a Riemannian sGS-type ADMM with a convergence guarantee is proposed.

Now, let us establish the equivalence between the original RADA-RGD (i.e., Algorithm 3) and its modified version (i.e., Algorithm 4).

The updates of x and y in Algorithm 3 are equivalent to (78) and (79) in Algorithm 4, respectively.

Proof Denote p⁡(x):=\argminp∈ℰ2​ℒ~1/βk​(x,p,yk). By [40, Lemma 1], we have

p⁡(x) =βkλ+βk​prox(λ+βk)​h∗​(𝒜⁡(x)+βk​yk)+λλ+βk​(𝒜⁡(x)+βk​yk),

which implies that

𝒜⁡(x)−p⁡(x)=βkλ+βk​(𝒜⁡(x)−λ​yk−prox(λ+βk)​h∗​(𝒜⁡(x)+βk​yk)). (80)

We first show the equivalence of the updates of y in the two algorithms. Note that p⁡(xk+1)=pk+1. Substituting (80) with x=xk+1 into (79) yields

yk+1= βkλ+βk​yk+1λ+βk​(𝒜⁡(xk+1)−prox(λ+βk)​h∗​(𝒜⁡(xk+1)+βk​yk)),

which is exactly the update (69), or equivalently, the update (14) in Algorithm 3.

We now show the equivalence of the updates of x in the two algorithms. Substituting (80) into ∇xℒ~1/βk​(x,p⁡(x),yk) yields

∇xℒ~1/βk​(x,p⁡(x),yk)= ∇f(x)+∇𝒜(x)⊤(yk+1βk(𝒜(x)−p(x)))

= ∇f(x)+1λ+βk∇𝒜(x)⊤(𝒜(x)+βkyk−prox(λ+βk)​h∗(𝒜(x)+βkyk)),

which, together with (18) and the Moreau decomposition (8), implies that

∇xℒ~1/βk​(x,p⁡(x),yk)=∇Φk​(x).

Therefore, noting that p⁡(xk,t)=pk,t, we can rewrite (78) as

xk,t+1=Rxk,t​(−ζk,t​grad​Φk​(xk,t)),

which is exactly the update of x in RADA-RGD.

6 Numerical Results

In this section, we report the numerical results of our proposed RADA algorithmic framework for solving SPCA, FPCA, and SSC problems. All the tests are implemented in MATLAB 2023b and evaluated on Apple M2 Pro CPU. Our code is available at https://github.com/XuMeng00124/RADAopt.

Table 2: Average performance comparison on SPCA.

RADA-RGD RALM-II ManPG-Ada

−Φ var iter cpu spar −Φ var iter cpu spar −Φ var iter cpu spar

μ Synthetic, d=1000,N=50,r=10

0.5 859.3 0.991 636 1.9 26.6 859.0 0.991 26 5.4 26.7 859.1 0.991 18366 11.1 26.5

0.75 810.2 0.979 2288 1.3 35.5 809.9 0.978 28 5.3 35.7 810.1 0.978 19581 12.0 35.5

1 765.8 0.965 494 1.5 42.1 765.5 0.964 30 5.2 42.3 765.4 0.964 24640 14.5 42.3

1.25 724.2 0.952 569 1.7 47.4 723.8 0.951 33 5.3 47.4 724.0 0.951 18271 12.0 47.4

1.5 685.8 0.936 527 1.6 51.8 685.7 0.936 36 5.4 51.9 685.3 0.936 13713 9.1 51.9

r Synthetic, d=2000,N=50,μ=3

8 1163.1 0.919 532 2.5 46.1 1161.6 0.923 61 16.2 47.6 1161.0 0.915 15004 12.4 49.2

10 1195.2 0.904 624 3.7 59.3 1193.6 0.896 62 21.1 59.8 1190.0 0.885 16399 19.2 61.0

12 1189.7 0.890 816 5.5 65.1 1188.7 0.891 70 29.6 64.3 1188.6 0.890 12614 22.6 65.2

14 1145.9 0.867 698 5.7 69.8 1144.5 0.869 64 31.5 69.0 1144.6 0.865 15179 46.8 69.9

16 1247.8 0.895 739 7.0 75.6 1245.0 0.892 50 28.5 76.1 1240.1 0.894 23799 98.3 76.0

μ coil-100, d=1024,N=100,r=10

1 531.0 0.960 268 1.1 47.2 530.8 0.959 34 5.6 47.4 530.4 0.959 17205 11.9 47.3

1.25 488.5 0.940 234 1.0 55.3 487.8 0.940 35 4.9 55.2 487.6 0.940 11865 8.9 55.2

1.5 453.5 0.923 266 1.2 61.3 451.1 0.917 39 5.1 61.7 448.5 0.911 11828 8.7 61.7

1.75 419.7 0.906 319 1.4 66.4 417.2 0.898 41 5.7 66.5 411.5 0.885 10475 7.7 66.5

2 388.2 0.890 344 1.5 70.6 380.5 0.866 43 5.6 70.2 373.7 0.848 10950 7.9 70.4

r coil-100, d=1024,N=100,μ=2

8 384.8 0.913 254 0.9 63.2 380.3 0.898 40 3.9 63.5 375.2 0.883 16313 7.1 64.2

10 388.2 0.890 344 1.5 70.6 380.5 0.866 43 5.7 70.2 373.7 0.848 10950 7.9 70.4

12 403.3 0.875 654 3.1 76.0 395.7 0.848 40 5.8 75.6 391.5 0.837 10463 10.7 75.7

14 399.9 0.860 636 3.5 80.5 391.2 0.829 41 8.1 79.3 384.5 0.814 14758 22.6 79.4

16 398.7 0.842 640 3.8 83.3 388.7 0.808 38 8.2 82.9 385.4 0.800 12004 25.3 82.8

Before reporting the results, we specify the default choices of some parameters in our tests. We choose ζmin=10−20,ζmax=1020, τ1=0.999,τ2=0.9,ρ=1.5, and λ=ε/(2​R) for RADA-PGD and RADA-RGD. We also set λ=ε/(2​R) for ARPGDA. For RADA-RGD, we set νk=2​Tk​R2​βk,c1=10−4, and η=0.1. The value of R can be easily computed for each specific problem and is given by μ​d​r, 1, and μ​N for problems (3), (4), and (6), respectively. The initial setting of the parameter β1 varies across different problem settings. For DSGM, ARPGDA, and RADMM, we apply the same line search strategy initialized with the BB stepsize as in RADA-RGD for updating x. Unless otherwise stated, we terminate RADA-PGD and RADA-RGD if they return an ε-RGS point. The stopping rules for the other algorithms will be specified later. We utilize the QR factorization and the standard projection as the retraction on the Stiefel manifold and the Grassmann manifold, respectively [1, 2, 6, 63]. The standard projection onto the Grassmann manifold is obtained by taking the m leading eigenvectors of the matrix to be projected. The Riemannian gradient on the Grassmann manifold is computed as in [63]. Finally, both the QR-based retraction on the Stiefel manifold and the projection-based retraction on the Grassmann manifold are globally defined for all tangent vectors. Therefore, Assumption 4.1 holds with ϱ=+∞ for both manifolds.

6.1 Results on SPCA

In this subsection, we compare our proposed RADA-RGD with RALM-II (i.e., LS-II in [83]) and ManPG-Ada [11], a more efficient version of ManPG, on synthetic datesets and the real dataset coil-100 [54], which contains 7200 image samples of 100 objects taken from different angles with d=1024. For RADA-RGD, we set β1=0.1​d​r and Tk≡10. The codes for RALM-II and ManPG-Ada are downloaded from https://github.com/miskcoo/almssn and https://github.com/chenshixiang/ManPG, respectively. For ManPG-Ada, we relax the line search parameter from 0.5 to 10−4 and adopt a nonmonotone line search strategy to improve stepsize selection and empirical performance. We terminate all three methods once an ε-ROS point is found. More specifically, we terminate RADA-RGD when

max{dist(0,gradf(xk)+projTxk​ℳ(∇𝒜(xk)⊤∂h∗(pk))),‖pk−𝒜(xk)‖}≤ε, (81)

where pk=yk+12 is defined in (20). For RALM-II, we apply the same criterion (81) but with pk=yk. Recall that ManPG solves the subproblem

Vk:=arg⁡minV∈Txk​ℳ​{⟨grad​f​(xk),V⟩+ℓf2​‖V‖2+h∗​(xk+V)}

at each iteration. Here, ℓf is the Lipschitz constant of ∇f. We terminate ManPG-Ada if max⁡{ℓf​‖Vk‖,‖Vk‖}≤ε. The maximum number of iterations is set as 50000. In each test, we generate the synthetic data matrix A following the method in [83] or randomly select 100 samples from the coil-100 dataset.

First, we set ε=10−4. The average results over 20 runs, each using either a different randomly generated synthetic matrix A or a different set of selected samples from the coil-100 dataset, are reported in Table 2. Here, Φ is defined in (21), “cpu” represents the cpu time in seconds, and “iter” denotes the outer iteration number. We also compare the normalized variance (denoted by “var”) and the sparsity of X (denoted by “spar”). The former is defined as ⟨A​A⊤,X^​X^⊤⟩/maxX∈𝒮⁡(d,r)⁡⟨A​A⊤,X​X⊤⟩, where X^ denotes the returned solution by the algorithm. The latter is defined as the percentage of entries with the absolute value less than 10−5.

From Table 2, we see that the solutions obtained by these algorithms are comparable in terms of the normalized variance and sparsity. Our proposed RADA-RGD always returns the best solution in terms of the value of Φ. Moreover, the proposed single-loop RADA-RGD algorithm is significantly more efficient than the compared nested-loop algorithms. The performance difference between RADA-RGD and RALM-II is consistent with the analysis in Section 5. Specifically, although RALM-II requires significantly fewer iterations compared to RADA-RGD, it requires a much higher level of accuracy than RADA-RGD when solving the subproblem. Consequently, RALM-II is less efficient than RADA-RGD for solving the SPCA problems considered.

Figure 1: Stationarity gap versus runtime for two representative instances of SPCA.

Next, we set ε=10−8 and tighten the subproblem stopping criterion of ManPG-Ada to 10−10 (otherwise, the stationarity gap stagnates and fails to reach the same accuracy level as RADA-RGD and RALM-II) to observe the stationarity gap versus the runtime of two representative instances in Figure 1. To ensure a fair and readable comparison, we record the runtime of RALM-II at each outer iteration and that of RADA-RGD and ManPG-Ada every 10 iterations. This is because RALM-II involves the fewest but most expensive outer iterations. As shown in Figure 1, RADA-RGD reduces the stationarity gap significantly faster than both RALM-II and ManPG-Ada, which is consistent with the results reported in Table 2.

6.2 Results on FPCA

For the FPCA problem (4), we generate the data points {ai}i=1N independently according to the standard Gaussian distribution as in [74] and divide them into m=N groups. We compare RADA-RGD with the single-loop algorithms ARPGDA [74], RADMM [40], and DSGM [4]. Note that the existing convergence guarantees for RADMM and DSGM do not apply to FPCA since the mapping 𝒜 is nonlinear. For RADA-RGD, we set β1=104​m2​r, Tk≡5, and ε=10−8. For ARPGDA, we set βk=103​m2​r/k1.5. For DSGM and RADMM, we set the smoothing parameter λk=10/k1/3. For RADMM, a significant effort is devoted to fine-tuning the penalty parameter in the augmented Lagrangian function, and the choice σk=10−7​k1.5 is found to work well in our tests. It should be noted that although this choice of σk violates the convergence condition in [40], the resulting RADMM demonstrates superior numerical performance. For the compared algorithms, a more relaxed criterion is employed: The algorithm is terminated if the improvement in Φ is less than 10−8 over 1000 consecutive iterations. The maximum number of iterations for all algorithms is set as 20000.

Table 3: Average performance comparison on FPCA.

RADA-RGD ARPGDA RADMM DSGM

−Φ cpu iter −Φ cpu iter −Φ cpu iter −Φ cpu iter

N d=1000,r=2

20 120.54 0.01 64 118.08 0.24 7715 120.18 0.20 6144 113.38 0.16 3123

25 98.88 0.01 77 97.00 0.24 6708 97.94 0.24 6885 97.01 0.19 3599

30 84.16 0.02 85 82.29 0.22 5764 83.11 0.30 7858 82.32 0.24 4276

35 73.69 0.02 102 72.16 0.33 7833 72.80 0.27 6289 69.74 0.24 3836

40 65.40 0.02 98 64.17 0.32 7257 64.77 0.25 5814 60.24 0.25 3824

d N=20,r=4

200 56.86 0.01 87 56.51 0.02 1579 56.59 0.02 1178 53.73 0.06 1890

400 104.24 0.01 73 103.29 0.07 3174 103.08 0.11 4771 103.39 0.11 2420

600 146.57 0.01 62 145.55 0.32 9975 146.20 0.22 6717 140.12 0.20 3821

800 193.26 0.01 53 191.28 0.38 9094 192.64 0.31 7248 192.25 0.32 5244

1000 237.18 0.02 75 235.59 0.51 10124 236.65 0.38 7423 224.79 0.38 5270

Table 3 presents the results on FPCA averaged over 20 runs with different randomly generated initial points. From this table, we observe that our proposed RADA-RGD consistently returns the best solutions in terms of the value of Φ while requiring the least computational time. Specifically, on average, our proposed RADA-RGD is more than 10 times faster than the other compared algorithms when d is large. We further plot the stationarity gap versus the number of gradient evaluations for two representative instances in Figure 2. As illustrated in the figure, RADA-RGD achieves a faster reduction in the stationarity gap than the other first-order methods. Moreover, the similarities and differences between RADA-RGD and RADMM, as discussed in Section 5.2, provide important insights into the superior performance of RADA-RGD. In particular, as shown in Algorithm 4, an additional update of p is performed before updaing x, which means that the Riemannian gradient used to update x in RADA-RGD involves the most up-to-date information of p. The superior performance of RADA-RGD over ARPGDA further demonstrates that performing multiple Riemannian gradient steps to ensure sufficient descent of the proposed value function Φk is much more efficient than simply performing a Riemannian gradient descent step on certain surrogate of F⁡(⋅,y) for some fixed y, as is done in ARPGDA and its Euclidean counterparts [76, 58, 53, 20].

Figure 2: Stationarity gap versus number of gradient evaluations for two representative instances of FPCA.

6.3 Results on SSC

In this subsection, we compare our proposed RADA-PGD and RADA-RGD with ARPGDA [74], DSGM [4], and NADMM [51] for solving the SSC problem (5). Note that NADMM is specifically designed for solving (5). The main difference between NADMM and RADMM is that NADMM performs a Euclidean projected gradient descent step to update x.

First, we perform tests on synthetic datasets, wherein the data points {ai}i=1N are independently generated according to the standard Gaussian distribution and Wi​j=|⟨ai,aj⟩| as in [52]. We take Q1=X1​X1⊤ as the initial point, where X1 consists of the m eigenvectors associated to the m smallest eigenvalues of L. For RADA, we set ε=10−4. For RADA-PGD, we set β1=N2​m, Tk≡1, and ζk,t=ℓk−1 with ℓk=(λ+βk)−1. For RADA-RGD, we set β1=N2​m and Tk≡3. For DSGM and NADMM, we set λk=10−2/k1/3. Additionally, we carefully tune the parameter βk in ARPGDA for each test to achieve its best possible performance. For NADMM, we set the penalty parameter σk in the augmented Lagrangian function as 0.1​k1/3, which enhances the efficiency of NADMM in our tests but does not satisfy the conditions of the theoretical guarantee in [51]. For the compared algorithms, a more relaxed criterion is employed: The algorithm is terminated if the improvement of the objective function Φ is less than 10−8 over 50 consecutive iterations.

Table 4: Average performance comparison on SSC (synthetic datasets).

RADA-PGD RADA-RGD ARPGDA NADMM DSGM

Φ cpu iter Φ cpu iter Φ cpu iter Φ cpu iter Φ cpu iter

m N=200,μ=0.005

2 1.976 0.1 72 1.976 0.2 67 1.980 0.6 580 1.978 0.6 838 1.983 8.0 4455

3 2.973 0.1 92 2.973 0.3 78 2.987 6.0 6494 2.977 0.8 1136 2.988 2.9 2168

4 3.970 0.2 104 3.970 0.2 77 3.988 5.6 6238 3.975 0.9 1324 3.989 4.5 3528

5 4.967 0.2 123 4.967 0.2 76 4.987 6.0 6446 4.973 1.2 1759 4.992 6.8 5375

6 5.965 0.2 141 5.965 0.3 86 5.987 5.8 6185 5.971 1.7 2421 5.994 9.8 7601

μ N=500,m=3

0.001 2.496 1.0 116 2.496 3.0 158 2.500 22.2 3757 2.498 3.8 1244 2.504 73.5 9198

0.002 2.985 1.1 110 2.985 2.0 93 2.990 46.5 7294 2.986 1.9 531 2.996 25.9 3084

0.005 3.007 1.7 152 3.007 3.2 143 3.010 3.3 503 3.007 3.0 798 3.009 29.7 3282

0.01 3.022 1.9 165 3.022 3.5 152 3.031 12.2 1796 3.022 13.9 3480 3.024 16.5 1681

0.02 3.052 2.6 227 3.052 4.1 164 3.057 3.2 436 3.052 17.1 3704 3.052 13.3 1134

Figure 3: Stationarity gap versus number of gradient evaluations for two representative instances of SSC.

Table 4 reports the results on SSC averaged over 20 runs with randomly generated data points, and Figure 3 illustrates the stationarity gap versus the number of gradient evaluations for two representative instances. From both the table and the figure, we observe that the proposed RADA-PGD and RADA-RGD consistently achieve the best values of Φ and effectively reduce the stationarity gap. Notably, RADA-PGD is the most efficient among all compared algorithms. Similar to RADA-RGD, RADA-PGD with Tk≡1 can also be interpreted as a Riemannian sGS-type ADMM. This interpretation sheds light on the superior performance of RADA-PGD over NADMM, given that sGS-type ADMM has demonstrated better performance than ADMM in many scenarios [10, 43, 44]. Furthermore, RADA-PGD demonstrates even better performance than RADA-RGD, especially when N=500. This can be attributed to two reasons. First, RADA-PGD only needs to compute the Euclidean gradient instead of the Riemannian gradient, which requires an additional projection onto the tangent space. Second, the parameter ℓk in (45) becomes (λ+βk)−1 for this specific problem, making it straightforward to compute. Consequently, we can use ℓk−1=λ+βk directly as an appropriate stepsize at the k-th iteration without any additional cost.

Next, to further evaluate the different algorithms, we compare their clustering performance on four real datasets (i.e., Wine, Iris, Glass, and Letter) from the UCI machine learning repository [35]. Due to the large size of the Letter dataset, we randomly select a subset of 1000 samples from the first 17 classes. The affinity matrix W is constructed using the Gaussian kernel as in [51], i.e., Wi​j=exp(−∥ai−aj∥2/κ). We fine-tune the parameters κ and μ for each dataset so that the solutions to the obtained optimization problems yield good clustering performance. We utilize the normalized mutual information (NMI) scores [66] to measure the clustering performance. Note that a higher NMI score indicates better clustering performance. The parameters and stopping criteria for all algorithms remain the same as those used in the tests on synthetic datasets, with the exception that ε=10−3 is used for both RADA-PGD and RADA-RGD. We also fine-tune the parameter βk for ARPGDA in each test.

Table 5: Performance comparison on SSC (UCI datasets).

RADA-PGD RADA-RGD ARPGDA DSGM NADMM

Wine: N=178,m=3,κ=1,μ=0.001

NMI 0.893 0.893 0.893 0.893 0.893

Φ 1.899 1.899 1.899 1.899 1.899

iter 148 115 932 1952 659

cpu 0.26 0.52 1.60 5.12 0.97

Iris: N=149,m=3,κ=0.2,μ=0.005

NMI 0.776 0.776 0.776 0.757 0.776

Φ 1.106 1.106 1.106 1.106 1.106

iter 161 115 320 1491 5995

cpu 0.21 0.44 0.59 3.17 4.82

Glass: N=213,m=3,κ=1,μ=0.0002

NMI 0.471 0.471 0.471 0.471 0.471

Φ 1.818 1.818 1.818 1.818 1.818

iter 120 63 268 8474 177

cpu 0.17 0.30 0.44 28.35 0.34

Letter: N=700,m=12,κ=1,μ=0.001

NMI 0.386 0.384 0.339 0.297 0.315

Φ 16.63 16.64 16.70 16.73 16.66

iter 271 222 4456 10000 10000

cpu 8.71 23.89 216.65 450.76 313.04

Table 5 reports the comparison of NMI scores, objective function values, and computational time for the four real datasets. From the table, we see that RADA-PGD and RADA-RGD consistently find solutions with the best quality in terms of NMI scores and objective function values for all datasets. For the Letter dataset, the proposed RADA-PGD enjoys a significantly lower computational time, while DSGM and NDAMM hit the preset maximum number of iterations of 10000.

7 Concluding Remarks

In this paper, we introduced a flexible RADA algorithmic framework based on a novel value function for solving a class of Riemannian NC-L minimax problems. Within this framework, we proposed two customized efficient algorithms called RADA-PGD and RADA-RGD, which perform one or multiple projected/Riemannian gradient descent steps and then a proximal gradient ascent step at each iteration. Theoretically, we proved that the proposed algorithmic framework can find both ε-RGS and ε-ROS points of the Riemannian NC-L minimax problem within 𝒪⁡(ε−3) iterations, achieving the best-known iteration complexity. We also showed how the proposed algorithmic framework can be adapted to solve an equivalent Riemannian nonsmooth reformulation of the minimax problem. As a by-product, we made the intriguing observation that RADA and RADA-RGD are closely related to, but essentially different from, RALM and RADMM, respectively. Computationally, we demonstrated the superiority of our proposed algorithms compared to existing state-of-the-art algorithms through extensive numerical results on SPCA, FPCA, and SSC. It is worth noting that our proposed RADA algorithmic framework and the associated algorithms are not derived from existing Euclidean methods; rather, the ideas behind their design even lead to new algorithms for NC-L minimax problems in the Euclidean setting.

Our work suggests several directions for future research. One direction is to design efficient algorithms with strong theoretical guarantees for solving more general Riemannian minimax problems. Another direction is to develop stochastic versions of our algorithms. It is also worthwhile to develop efficient algorithms for broader classes of Riemannian nonsmooth composite problems by suitably exploiting their intrinsic minimax structure.

The authors are grateful to the Area Editor, the Associate Editor, and the anonymous reviewers for their valuable comments and suggestions, which have helped improve the quality of the manuscript.

References

[1] P.-A. Absil, R. Mahony, and R. Sepulchre (2008) Optimization algorithms on matrix manifolds. Princeton Univ. Press. Cited by: §2, §6.

[2] P.-A. Absil and J. Malick (2012) Projection-like retractions on matrix manifolds. SIAM J. Optim. 22 (1), pp. 135–158. Cited by: §6.

[3] J. Barzilai and J. M. Borwein (1988) Two-point step size gradient methods. IMA J. Numer. Anal. 8 (1), pp. 141–148. Cited by: §4.

[4] A. Beck and I. Rosset (2023) A dynamic smoothing technique for a class of nonsmooth optimization problems on manifolds. SIAM J. Optim. 33 (3), pp. 1473–1493. Cited by: item ∙ , item ∙ , §1.2, Table 1, §3.2, §3.2, §5.1, §6.2, §6.3.

[5] A. Beck (2017) First-Order Methods in Optimization. SIAM. Cited by: §1.2, Theorem 2.1, §2.

[6] T. Bendokat, R. Zimmermann, and P.-A. Absil (2024) A Grassmann manifold handbook: basic geometry and computational aspects. Adv. Comput. Math. 50 (1), pp. 1–51. Cited by: item 3., §6.

[7] P. B. Borckmans, S. E. Selvan, N. Boumal, and P.-A. Absil (2014) A Riemannian subgradient algorithm for economic dispatch with valve-point effect. J. Comput. Appl. Math. 255, pp. 848–866. Cited by: §1.2.

[8] N. Boumal, P.-A. Absil, and C. Cartis (2019) Global rates of convergence for nonconvex optimization on manifolds. IMA J. Numer. Anal. 39 (1), pp. 1–33. Cited by: §4.1, §4.1, Proof 4.3.

[9] N. Boumal (2023) An introduction to optimization on smooth manifolds. Cambridge Univ. Press. Cited by: §2.

[10] L. Chen, D. Sun, and K.-C. Toh (2017) An efficient inexact symmetric Gauss–Seidel based majorized ADMM for high-dimensional convex composite conic programming. Math. Program. 161, pp. 237–270. Cited by: §5.2, §6.3.

[11] S. Chen, S. Ma, A. M.-C. So, and T. Zhang (2020) Proximal gradient method for nonsmooth optimization over the Stiefel manifold. SIAM J. Optim. 30 (1), pp. 210–239. Cited by: item ∙ , §1.2, Table 1, §1, §4.1, §6.1.

[12] S. Chen, S. Ma, A. M.-C. So, and T. Zhang (2024) Nonsmooth optimization over the Stiefel manifold and beyond: proximal gradient method and recent variants. SIAM Rev. 66 (2), pp. 319–352. Cited by: §1.2, §4.1.

[13] K. Deng, J. Hu, J. Wu, and Z. Wen (2025) Oracle complexities of augmented Lagrangian methods for nonsmooth manifold optimization. Math. Oper. Res.. Note: doi: 10.1287/moor.2024.0498 Cited by: item ∙ , item ∙ , §1.2, Table 1, §3.2, §5.1, §5.1, §5.

[14] K. Deng and Z. Peng (2023) A manifold inexact augmented Lagrangian method for nonsmooth optimization on Riemannian submanifolds in Euclidean space. IMA J. Numer. Anal. 43 (3), pp. 1653–1684. Cited by: item ∙ , §1.2, §5.1, §5.1, §5.

[15] D. Drusvyatskiy and C. Paquette (2019) Efficiency of minimizing compositions of convex functions and smooth maps. Math. Program. 178, pp. 503–558. Cited by: §4.1.

[16] B. Gao, X. Liu, X. Chen, and Y. Yuan (2018) A new first-order algorithmic framework for optimization problems with orthogonality constraints. SIAM J. Optim. 28 (1), pp. 302–332. Cited by: §4.

[17] D. Hajinezhad and M. Hong (2019) Perturbed proximal primal-dual algorithm for nonconvex nonsmooth optimization. Math. Program. 176 (1), pp. 207–245. Cited by: §5.1.

[18] A. Han, B. Mishra, P. Jawanpuria, and J. Gao (2023) Nonconvex-nonconcave min-max optimization on Riemannian manifolds. Trans. Mach. Learn. Res., pp. . Cited by: §1.2, §1.2.

[19] A. Han, B. Mishra, P. Jawanpuria, P. Kumar, and J. Gao (2023) Riemannian Hamiltonian methods for min-max optimization on manifolds. SIAM J. Optim. 33 (3), pp. 1797–1827. Cited by: §1.2, §1.2.

[20] J. He, H. Zhang, and Z. Xu (2024) An approximation proximal gradient algorithm for nonconvex-linear minimax problems with nonconvex nonsmooth terms. J. Glob. Optim. 90, pp. 73–92. Cited by: item ∙ , item ∙ , §1.2, §3.1, §3.1, §3.1, §3.2, §3.2, §4, §6.2.

[21] S. Hosseini, W. Huang, and R. Yousefpour (2018) Line search algorithms for locally Lipschitz functions on Riemannian manifolds. SIAM J. Optim. 28 (1), pp. 596–619. Cited by: §1.2.

[22] S. Hosseini and A. Uschmajew (2017) A Riemannian gradient sampling algorithm for nonsmooth optimization on manifolds. SIAM J. Optim. 27 (1), pp. 173–189. Cited by: §1.2.

[23] J. Hu, X. Liu, Z. Wen, and Y. Yuan (2020) A brief introduction to manifold optimization. J. Oper. Res. Soc. China 8, pp. 199–248. Cited by: §1.

[24] J. Hu, A. Milzarek, Z. Wen, and Y. Yuan (2018) Adaptive quadratically regularized Newton method for Riemannian optimization. SIAM Journal on Matrix Analysis and Applications 39 (3), pp. 1181–1207. Cited by: §4.

[25] X. Hu, N. Xiao, X. Liu, and K. Toh (2024) A constraint dissolving approach for nonsmooth optimization over the Stiefel manifold. IMA J. Numer. Anal. 44 (6), pp. 3717–3748. Cited by: §1.2.

[26] F. Huang and S. Gao (2023) Gradient descent ascent for minimax problems on Riemannian manifolds. IEEE Trans. Pattern Anal. Mach. Intell. 45 (7), pp. 8466–8476. Cited by: §1.2, §1.2.

[27] W. Huang and K. Wei (2022) Riemannian proximal gradient methods. Math. Program. 194 (1), pp. 371–413. Cited by: §1.2, §1, §4.1.

[28] W. Huang and K. Wei (2023) An inexact Riemannian proximal gradient method. Comput. Optim. Appl. 85 (1), pp. 1–32. Cited by: §1.2, §4.1.

[29] B. Iannazzo and M. Porcelli (2018) The Riemannian Barzilai-Borwein method with nonmonotone line search and the matrix geometric mean computation. IMA J. Numer. Anal. 38 (1), pp. 495–517. Cited by: §4.

[30] B. Jiang and Y.-H. Dai (2015) A framework of constraint preserving update schemes for optimization on Stiefel manifold. Math. Program. 153 (2), pp. 535–575. Cited by: §4.

[31] B. Jiang and Y.-F. Liu (2023) A Riemannian exponential augmented Lagrangian method for computing the projection robust Wasserstein distance. In Proc. Adv. Neural Inf. Process. Syst., Vol. 36, pp. 79999–80023. Cited by: §4.

[32] C. Jin, P. Netrapalli, and M. Jordan (2020) What is local optimality in nonconvex-nonconcave minimax optimization?. In Proc. Int. Conf. Mach. Learn., pp. 4880–4889. Cited by: §1.2.

[33] I. T. Jolliffe, N. T. Trendafilov, and M. Uddin (2003) A modified principal component technique based on the LASSO. J. Comput. Graph. Stat. 12 (3), pp. 531–547. Cited by: item 1..

[34] M. Jordan, T. Lin, and E.-V. Vlatakis-Gkaragkounis (2022) First-order algorithms for min-max optimization in geodesic metric spaces. In Proc. Adv. Neural Inf. Process. Syst., Vol. 35, pp. 6557–6574. Cited by: §1.2, §1.2.

[35] M. Kelly, R. Longjohn, and K. Nottingham (2007) The UCI Machine Learning Repository. Note: https://archive.ics.uci.edu Cited by: §6.3.

[36] W. Kong and R. D. C. Monteiro (2021) An accelerated inexact proximal point method for solving nonconvex-concave min-max problems. SIAM J. Optim. 31 (4), pp. 2558–2585. Cited by: §1.2.

[37] W. Kong and R. D. C. Monteiro (2024) Global complexity bound of a proximal ADMM for linearly constrained nonseparable nonconvex composite programming. SIAM J. Optim. 34 (1), pp. 201–224. Cited by: §5.1.

[38] A. Kovnatsky, K. Glashoff, and M. M. Bronstein (2016) MADMM: A generic algorithm for non-smooth optimization on manifolds. In Proc. Comput. Vis. ECCV, pp. 680–696. Cited by: §1.2.

[39] R. Lai and S. Osher (2014) A splitting method for orthogonality constrained problems. J. Sci. Comput. 58, pp. 431–449. Cited by: §1.2.

[40] J. Li, S. Ma, and T. Srivastava (2024) A Riemannian alternating direction method of multipliers. Math. Oper. Res. 50 (4), pp. 3222–3242. Cited by: item ∙ , item ∙ , §1.2, Table 1, §1, §3.2, §5.2, §5.2, §5.2, Proof 5.2, §5, §6.2.

[41] J. Li, L. Zhu, and A. M.-C. So (2025) Nonsmooth nonconvex-nonconcave minimax optimization: primal-dual balancing and iteration complexity analysis. Math. Program. 214, pp. 591–641. Cited by: §1.2, §3.2.

[42] X. Li, S. Chen, Z. Deng, Q. Qu, Z. Zhu, and A. M.-C. So (2021) Weakly convex optimization over Stiefel manifold using Riemannian subgradient-type methods. SIAM J. Optim. 31 (3), pp. 1605–1634. Cited by: §1.2, Table 1, §1.

[43] X. Li, D. Sun, and K.-C. Toh (2016) A Schur complement based semi-proximal ADMM for convex quadratic conic programming and extensions. Math. Program. 155 (1), pp. 333–373. Cited by: §5.2, §6.3.

[44] X. Li, D. Sun, and K.-C. Toh (2019) A block symmetric Gauss-Seidel decomposition theorem for convex composite quadratic programming and its applications. Math. Program. 175, pp. 395–418. Cited by: §5.2, §6.3.

[45] T. Lin, C. Jin, and M. I. Jordan (2020) Near-optimal algorithms for minimax optimization. In Proc. Conf. Learn. Theory, pp. 2738–2779. Cited by: §1.2.

[46] T. Lin, C. Jin, and M. Jordan (2020) On gradient descent ascent for nonconvex-concave minimax problems. In Proc. Int. Conf. Mach. Learn., pp. 6083–6093. Cited by: item ∙ , §1.2, §3.2.

[47] C. Liu and N. Boumal (2020) Simple algorithms for optimization on Riemannian manifolds with constraints. Appl. Math. Optim. 82 (3), pp. 949–981. Cited by: §1.2.

[48] H. Liu, A. M.-C. So, and W. Wu (2019) Quadratic optimization with orthogonality constraint: explicit Łojasiewicz exponent and linear convergence of retraction-based line-search and stochastic variance-reduced gradient methods. Math. Program. 178 (1), pp. 215–262. Cited by: §4.1.

[49] X. Liu, N. Xiao, and Y. Yuan (2024) A penalty-free infeasible approach for a class of nonsmooth optimization problems over the Stiefel manifold. J. Sci. Comput. 99 (2), pp. 30. Cited by: §1.2.

[50] Y.-F. Liu, T.-H. Chang, M. Hong, Z. Wu, A. M.-C. So, E. A. Jorswieck, and W. Yu (2024) A survey of recent advances in optimization methods for wireless communications. IEEE J. Sel. Areas Commun. 42 (11), pp. 2992–3031. Cited by: §1.

[51] C. Lu, J. Feng, Z. Lin, and S. Yan (2018) Nonconvex sparse spectral clustering by alternating direction method of multipliers and its convergence analysis. In Proc. AAAI Conf. Artif. Intell., Vol. 32. Cited by: item 3., item ∙ , §6.3, §6.3, §6.3.

[52] C. Lu, S. Yan, and Z. Lin (2016) Convex sparse spectral clustering: Single-view to multi-view. IEEE Trans. Image Process. 25 (6), pp. 2833–2843. Cited by: item 3., §6.3.

[53] S. Lu, I. Tsaknakis, M. Hong, and Y. Chen (2020) Hybrid block successive approximation for one-sided non-convex min-max problems: algorithms and applications. IEEE Trans. Signal Process. 68, pp. 3676–3691. Cited by: item ∙ , §1.2, Table 1, §3.1, §3.1, §3.1, §4, §6.2.

[54] S. A. Nene, S. K. Nayar, and H. Murase (1996) Columbia object image library (COIL-100). Technical Report CUCS-006-96. Cited by: §6.1.

[55] Y. Nesterov (2018) Lectures on convex optimization. Vol. 137, Springer. Cited by: Proof 4.3.

[56] M. Nouiehed, M. Sanjabi, T. Huang, J. D. Lee, and M. Razaviyayn (2019) Solving a class of non-convex min-max games using iterative first order methods. In Proc. Adv. Neural Inf. Process. Syst., Vol. 32, pp. 14934–14942. Cited by: §1.2.

[57] D. M. Ostrovskii, A. Lowy, and M. Razaviyayn (2021) Efficient search of first-order Nash equilibria in nonconvex-concave smooth min-max problems. SIAM J. Optim. 31 (4), pp. 2508–2538. Cited by: §1.2.

[58] W. Pan, J. Shen, and Z. Xu (2021) An efficient algorithm for nonconvex-linear minimax optimization problem and its application in solving weighted maximin dispersion problem. Comput. Optim. Appl. 78, pp. 287–306. Cited by: item ∙ , item ∙ , §1.2, Table 1, §3.1, §3.1, §3.1, §3.2, §3.2, §4, §6.2.

[59] Z. Peng, W. Wu, J. Hu, and K. Deng (2023) Riemannian smoothing gradient type algorithms for nonsmooth optimization problem on compact Riemannian submanifold embedded in Euclidean space. Appl. Math. Optim. 88 (3), pp. 85. Cited by: §1.2, §5.1.

[60] H. Rafique, M. Liu, Q. Lin, and T. Yang (2022) Weakly-convex-concave min-max optimization: provable algorithms and applications in machine learning. Optim. Methods Softw. 37 (3), pp. 1087–1121. Cited by: §1.2.

[61] R. T. Rockafellar and R. J.-B. Wets (2009) Variational analysis. Vol. 317, Springer Sci. Business Media. Cited by: §2.

[62] S. Samadi, U. Tantipongpipat, J. H. Morgenstern, M. Singh, and S. Vempala (2018) The price of fair PCA: one extra dimension. In Proc. Adv. Neural Inf. Process. Syst., Vol. 31, pp. 10999–11010. Cited by: item 2..

[63] H. Sato and T. Iwai (2014) Optimization algorithms on the Grassmann manifold with application to matrix eigenvalue problems. Jpn. J. Ind. Appl. Math. 31, pp. 355–400. Cited by: item 3., §6.

[64] J. Shen, A. J. Davis, D. Lu, and Z. Bai (2025) Hidden convexity of fair PCA and fast solver via eigenvalue optimization. arXiv prepr. arXiv:2503.00299. Cited by: item 2..

[65] J. Shen, Z. Wang, and Z. Xu (2023) Zeroth-order single-loop algorithms for nonconvex-linear minimax problems. J. Glob. Optim. 87 (2), pp. 551–580. Cited by: §1.2, §3.2.

[66] A. Strehl and J. Ghosh (2002) Cluster ensembles–A knowledge reuse framework for combining multiple partitions. J. Mach. Learn. Res. 3 (), pp. 583–617. Cited by: §6.3.

[67] K. Sun and X. A. Sun (2024) Dual descent augmented Lagrangian method and alternating direction method of multipliers. SIAM J. Optim. 34 (2), pp. 1679–1707. Cited by: §3.1.

[68] K. K. Thekumparampil, P. Jain, P. Netrapalli, and S. Oh (2019) Efficient algorithms for smooth minimax optimization. In Proc. Adv. Neural Inf. Process. Syst., Vol. 32, pp. 12659–12670. Cited by: §1.2.

[69] L. Tian and A. M.-C. So (2024) No dimension-free deterministic algorithm computes approximate stationarities of Lipschitzians. Math. Program. 208 (1), pp. 51–74. Cited by: §3.2.

[70] Z. Wang, B. Liu, S. Chen, S. Ma, L. Xue, and H. Zhao (2022) A manifold proximal linear method for sparse spectral clustering with application to single-cell RNA sequencing data analysis. INFORMS J. Optim. 4 (2), pp. 200–214. Cited by: item 3., §1.2, Table 1.

[71] Z. Wen and W. Yin (2013) A feasible method for optimization with orthogonality constraints. Math. Program. 142 (1), pp. 397–434. Cited by: §4.

[72] Z. Wu, B. Jiang, Y.-F. Liu, M. Shao, and Y.-H. Dai (2024) Efficient CI-based one-bit precoding for multiuser downlink massive MIMO systems with PSK modulation. IEEE Trans. Wireless Commun. 23 (5), pp. 4861–4875. Cited by: §1.2.

[73] N. Xiao, X. Liu, and Y. Yuan (2021) Exact penalty function for ℓ2,1 norm minimization over the Stiefel manifold. SIAM J. Optim. 31 (4), pp. 3097–3126. Cited by: §1.

[74] M. Xu, B. Jiang, W. Pu, Y.-F. Liu, and A. M.-C. So (2024) An efficient alternating Riemannian/projected gradient descent ascent algorithm for fair principal component analysis. In Proc. IEEE Int. Conf. Acoust., Speech, Signal Process., Vol. , pp. 7195–7199. Cited by: item 2., item ∙ , §1.2, Table 1, §3.1, §3.2, §3.2, §4, §6.2, §6.3.

[75] Z. Xu, Z. Wang, J. Shen, and Y. Dai (2024) Derivative-free alternating projection algorithms for general nonconvex-concave minimax problems. SIAM J. Optim. 34 (2), pp. 1879–1908. Cited by: §1.2.

[76] Z. Xu, H. Zhang, Y. Xu, and G. Lan (2023) A unified single-loop alternating gradient projection algorithm for nonconvex-concave and convex-nonconcave minimax problems. Math. Program. 201 (1), pp. 635–706. Cited by: item ∙ , §1.2, §3.1, §3.1, §3.1, §3.2, §3.2, §4, §6.2.

[77] J. Yang, A. Orvieto, A. Lucchi, and N. He (2022) Faster single-loop algorithms for minimax optimization without strong concavity. In Proc. Int. Conf. Artif. Intell. Stat., pp. 5485–5517. Cited by: item ∙ , §1.2, §3.2.

[78] J. Yang, S. Zhang, N. Kiyavash, and N. He (2020) A Catalyst framework for minimax optimization. In Proc. Adv. Neural Inf. Process. Syst., Vol. 33, pp. 5667–5678. Cited by: §1.2.

[79] G. Zalcberg and A. Wiesel (2021) Fair principal component analysis and filter design. IEEE Trans. Signal Process. 69, pp. 4835–4842. Cited by: item 2..

[80] C. Zhang, X. Chen, and S. Ma (2024) A Riemannian smoothing steepest descent method for non-Lipschitz optimization on embedded submanifolds of ℝn. Math. Oper. Res. 49 (3), pp. 1710–1733. Cited by: §1.2, §5.1.

[81] J. Zhang, P. Xiao, R. Sun, and Z. Luo (2020) A single-loop smoothed gradient descent-ascent algorithm for nonconvex-concave min-max problems. In Proc. Adv. Neural Inf. Process. Syst., Vol. 33, pp. 7377–7389. Cited by: §1.2.

[82] P. Zhang, J. Zhang, and S. Sra (2023) Sion’s minimax theorem in geodesic metric spaces and a Riemannian extragradient algorithm. SIAM J. Optim. 33 (4), pp. 2885–2908. Cited by: §1.2, §1.2.

[83] Y. Zhou, C. Bao, C. Ding, and J. Zhu (2023) A semismooth Newton based augmented Lagrangian method for nonsmooth optimization on matrix manifolds. Math. Program. 201 (1), pp. 1–61. Cited by: item ∙ , item ∙ , §1.2, §5.1, §5.1, §5, §6.1, §6.1.

[84] H. Zou and L. Xue (2018) A selective overview of sparse principal component analysis. Proc. IEEE 106 (8), pp. 1311–1320. Cited by: item 1..