Escaping Saddle Point Efficiently in Minimax and Bilevel Optimizations

1 Introduction

Hierarchical optimization (including minimax and bilevel optimization) is a popular and important optimization framework which is applied to a wide range of machine learning problems, such as Generative Adversarial Net [Goodfellow et al., 2014], adversarial training [Madry et al., 2018], multi-agent reinforcement learning [Wai et al., 2018], metalearning [Franceschi et al., 2018; Bertinetto et al., 2018] and hyperparameter optimization [Shaban et al., 2019; Feurer and Hutter, 2019]. In this paper, we study the following stochastic hierarchical optimization problem

where the upper-level function f(x, y∗ (x)) is smooth and possibly nonconvex, and the lower-level function g(x, y) is smooth and strongly-convex w.r.t. variable y so that y∗ (x) and Φ(x) can be well defined. ξ and ζ are samples drawn from data distribution D and D′ . Stochastic problem is a general form that covers a couple of optimization tasks, including online optimization and finite-sum optimization. When g(x, y) = -f(x, y), the above bilevel optimization problem is reduced to a standard minimax optimization which can be rewritten as Eq. (2)

where Y is a convex domain (not required to be compact). The loss function f(x, y) is smooth, nonconvex w.r.t. x and strongly-concave w.r.t. y.

1.1 Minimax Optimzation

Recently, there are plenty of works studying minimax optimization problem in a variety of research fields in machine learning. Many deterministic and stochastic algorithms with asymptotic or non-asymptotic convergence analysis have been developed, such as Gradient Descent Ascent (GDA) [Du and Hu, 2019; Nemirovski, 2004] and Stochastic Gradient Descent Ascent (SGDA) [Lin et al., 2020b]. Some algorithms adopt a single loop structure [Heusel et al., 2017; Lin et al., 2020b; Xu et al., 2023] while the others use a nested loop to update y more frequently so that they can obtain a better estimation of the maximum y∗ (x) [Jin et al., 2020; Nouiehed et al., 2019].

Besides, some algorithms have been proposed to improve the theoretical results of minimax optimization, such as SREDA [Luo et al., 2020] and Acc-MDA [Huang et al., 2022] which take advantage of variance reduction to accelerate the convergence rate and reduce the gradient complexity. Moreover, on deterministic setting some recently proposed algorithms [Lin et al., 2020a] have already matched the optimal lower bound [Zhang et al., 2021].

However, most of these works only consider the criterion of finding first-order stationary point. In nonconvex setting, convergence to first-order stationary point is not always satisfactory because a first-order stationary point could be a local minimum, saddle point or even local maximum. Therefore, second-order stationary point that reaches local minimum becomes a popular and important issue in nonconvex optimization. Since finding global minimum in nonconvex optimization is usually an NP-hard problem [Hillar and Lim, 2013], in some situations we attempt to find a local minimum instead. Moreover, in some machine learning tasks such as tensor decomposition [Ge et al., 2015], matrix sensing [Bhojanapalli et al., 2016; Park et al., 2017], and matrix completion [Ge et al., 2016], finding local minimum is equivalent to finding global minimum, which makes second-order stationary point more crucial.

Therefore, we are motivated to study the second-order optimality for minimax optimization which captures the local minimum of Φ(x). We will discuss the relation between the second-order stationary point of Φ(x) and the local equilibrium of f(x, y) in Section 3.2. In Section 3.1 we can see that under certain conditions the primal function Φ(x) is twice differentiable. An O(ϵ, ϵH) second-order stationary point satisfies ∥∇Φ(x)∥ ≤ O(ϵ) and λmin(∇2Φ(x)) ≥ -ϵH where λmin(·) means the smallest eigenvalue.

Although several recent works have been proposed to study the second-order stationary point for nonconvex-stronglyconcave minimax optimization based on cubic-regularized gradient descent ascent [Chen et al., 2021; Luo and Chen, 2022] or perturbed gradient [Huang et al., 2025], they are only adaptive to deterministic gradient oracle and finite-sum problem. The study of the second-order stationary point for stochastic nonconvex minimax problem where the full gradient is not available is still limited. A comparison between related minimax algorithms is demonstrated in Table 1.

Thus, to fill in this gap, we propose a new algorithm named Perturbed Recursive Gradient Descent Ascent (PRGDA) to search second-order stationary point for stochastic nonconvex problem (2). To our best knowledge, PRGDA is the first algorithm that is guaranteed to obtain second-order stationary point for stochastic nonconvex minimax optimization problems. Furthermore, our method is a pure first-order algorithm that only requires the computation of gradient oracle, which makes our method more efficient to implement. We also provide the analysis to show that the gradient complexity of our algorithm is ˜O(κ3 ϵ-3 ) to achieve O(ϵ, √ρΦϵ) second-order stationary point where κ and ρΦ are defined in Section 3.1, which matches the best result of finding first-order stationary point for the same stochastic nonconvex minimax problem.

1.2 Bilevel Optimization

Recently, many bilevel algorithms are proposed, such as deterministic algorithms AID-BiO and ITD-BiO [Ji et al., 2021], and stochastic algorithms BSA [Ghadimi and Wang, 2018] and StocBiO [Ji et al., 2021]. These methods are proposed to improve the convergence analysis of bilevel optimization since most earlier works [Domke, 2012; Pedregosa, 2016] only provide the asymptotic convergence analysis.

StocBiO algorithm [Ji et al., 2021] is a recent work to solve stochastic nonconvex-strongly-convex bilevel optimization via AID. In this paper, we also study the convergence of our method under this condition where Φ(x) is stochastic and probably nonconvex. According to previous studies of bilevel optimization, when f(x, y) and g(x, y) are differentiable and g(x, y) is strongly-convex with respect to y, Φ(x) is also differentiable and automatically ∥∇Φ(x)∥ ≤ ϵ is a criterion of first-order stationary point. Notice that in [Ji et al., 2021] ∥∇Φ(x)∥2 ≤ ϵ is used as the criterion. In this paper, we will uniformly adopt ∥∇Φ(x)∥ ≤ ϵ as the convergence criterion. More recently, many stochastic algorithms with variance reduction are proposed, such as RSVRB [Guo et al., 2021], SUSTAIN [Khanduri et al., 2021], MRBO and VRBO [Yang et al., 2021]. The gradient complexity of bilevel optimization is enhanced to O(ϵ-3 ), which is the best theoretical result as far as we know. StocBiO with iNEON [Huang et al., 2025] is a recent work that combines StocBiO with pure first-order method inexact negative curvature originated from noise (iNEON) to escape saddle point and find second-order stationary point for nonconvex-strongly-convex bilevel optimization.

Although these works are proposed to improve the performance of algorithms for bilevel optimization, the complexity of current methods that achieve second-order stationary point are still high. Actually, the complexity of StocBiO with iNEON is even higher than the standard StocBiO algorithm in order to find a local minimum with high probability. Thus, to fill these gap, we are motivated to propose an accelerated algorithm with variance reduction that requires lower complexity to find second-order stationary point for stochastic nonconvex-strongly-convex bilevel optimization.

The comparison between our method and related works to find O(ϵ) first-order stationary point or O(ϵ, √ρΦϵ) secondorder stationary point is shown in Table 2, where Gc(f, ϵ) and Gc(g, ϵ) are the numbers of gradient evaluations of function f(x, y) and g(x, y) respectively. The last column represents whether the algorithm can escape saddle point and find local minimum. Notation ˜O(·) hides the logarithm term. StocBiO with iNEON and our PRGDA algorithm involve a logarithm term in the complexity because they converge to second-order stationary point with high probability.

From Table 2 we can see our PRGDA algorithm improves the gradient complexity Gc(f, ϵ) and Gc(g, ϵ) of StocBiO with iNEON algorithm significantly and matches state-ofthe-art complexity O(ϵ-3 ), which is one of the most important contribution of this paper.

1.3 Contributions

We summarize our main contributions as follows:

• We propose a new PRGDA algorithm which is the first algorithm to reach second-order stationary point for stochastic nonconvex minimax optimization problem. Our method is pure first-order and does not require any calculation of second-order derivatives. Our method does not involve nested loops either, which makes it more efficient to implement.
• We prove that the gradient complexity of our algorithm is ˜O(κ3 ϵ-3 ) to achieve O(ϵ, √ϵ) second-order stationary point in stochastic nonconvex minimax optimization, which matches the best result of finding first-order stationary point in the same problem.
• PRGDA can also be applied to nonconvex bilevel optimization and we can prove that the gradient complexity is Gc(f, ϵ) = ˜O(κ3 ϵ-3 ) and Gc(g, ϵ) = ˜O(κ7 ϵ-3 ) to find O(ϵ, √ϵ) second-order stationary point in stochastic nonconvex bilevel optimization, which outperforms the previous best theoretical results and matches stateof-the-art to find first-order stationary point.

2 RelatedWork

2.1 Stochastic Minimax Optimization

Many algorithms are proposed to solve stochastic nonconvexstrongly-concave minimax problem, including intuitive methods SGDmax [Jin et al., 2020] and Stochastic Gradient Descent Ascent (SGDA) [Lin et al., 2020b]. Recently, some methods integrate variance reduction with minimax problem to accelerate the convergence, such as Stochastic Recursive gradiEnt Descent Ascent (SREDA) [Luo et al., 2020], Hybrid Variance-Reduced SGD [Tran Dinh et al., 2020] and Acc-MDA [Huang et al., 2022]. There are also some works that study the weakly-convex concave minimax optimization such as [Rafique et al., 2022] and [Yan et al., 2020]. More related to this work, Cubic-Regularized Gradient Descent Ascent (Cubic-GDA) [Chen et al., 2021] and Minimax Cubic Newton (MCN) [Luo and Chen, 2022] are two recent algorithms that can reach the second-order stationary point in nonconvex-strongly-concave minimax optimization.

2.2 Perturbed Gradient Descent

Perturbed Gradient Descent (PGD) [Jin et al., 2017] was proposed to find second-order stationary point for nonconvex optimization which introduces a perturbation under specific condition. It is a deterministic gradient based algorithm and only involves first-order oracle. Perturbed Gradient Descent algorithm consists of two phases, a descent phase and an escaping phase. In the descent phase, the algorithm runs gradient descent to make the function value decrease until the magnitude of the gradient is smaller than a certain threshold. In the escaping phase, it first introduces a perturbation drawn from a uniform distribution on the ball B0(r) with center 0 and radius r. After certain iterations of gradient descent, if the function value is reduced by a significant threshold then it indicates that the algorithm escapes a saddle point and it will do the descent phase again. Otherwise, it can be proven that the point where the last descent terminates is secondorder stationary with high probability. To extend PGD to the stochastic setting and incorporate it with variance reduction, SSRGD [Li, 2019] was proposed to reach second-order stationary point with stochastic first-order oracle (SFO) of O(ϵ-3.5 ). After that Pullback algorithm [Chen et al., 2022] was proposed to improve the complexity to O(ϵ-3 ).

2.3 Cubic-GDA andMinimax Cubic Newton

Cubic-Regularized Gradient Descent Ascent (Cubic-GDA) [Chen et al., 2021] and Minimax Cubic Newton (MCN) [Luo and Chen, 2022] are two recent algorithms that can reach the second-order stationary point in nonconvex-strongly-concave minimax optimization. Both of these two algorithms are inspired by cubic regularization and designed for deterministic problem. Cubic regularization was first proposed in [Nesterov and Polyak, 2006] which is a standard method that converges to second-order stationary point in conventional nonconvex optimization. However, Cubic-GDA and MCN are both designed for deterministic problem and neither of them works for the stochastic minimax problem (2) considered in this paper. Therefore, we are motivated to propose an algorithm that is suitable for the stochastic problem. Besides, Cubic-GDA and MCN involves the calculation of secondorder oracle or Hessian vector product while our method only requires the first-order information, which indicates that our method is more efficient to implement because the computation cost of Hessian matrix could be high.

2.4 StocBiOwith iNEON

In [Huang et al., 2025], algorithms for both minimax and bilevel optimization are proposed to find second-order stationary point. However, for minimax optimization only the deterministic problem is studied. For bilevel optimization, the stochastic problem is considered and the StocBiO with

Algorithm 1 Perturbed Recursive Gradient Descent Ascent

Input: initial value x0, y0

Parameter: stepsize η and ηH, perturbation radius r, escaping phase threshold tthres, average movement ¯D, tolerance ϵ, maximum iteration T .

1: Set escape = false, s = 0, esc = 0. 2: for t = 0, 1, . . . , T -1 do 3: Minimax: Update yt+1, vt, ut fromAlgorithm 2. 4: Bilevel: Update yt+1, vt, ut fromAlgorithm 3. 5: if escape = false then 6: if ∥vt∥ ≥ ϵ then 7: Update xt+1 = xt - (η/∥vt∥)vt. 8: else 9: Let ms = t, s = s + 1. 10: Set escape = true, esc = 0. 11: Draw perturbation ξ ∼ B0(r). 12: update xt+1 = xt + ξ. 13: end if 14: else 15: Compute D = ∑tj=ms+1 η2H ∥vj∥2 . 16: if D > (t - ms) ¯D then 17: Set ηt s.t. ∑tj=ms+1 η2 j ∥vj∥2 = (t - ms) ¯D. 18: Update xt+1 = xt - ηtvt. Set escape = false. 19: else 20: Set ηt = ηH. 21: Update xt+1 = xt - ηtvt, esc = esc + 1. 22: Return xms if esc = tthres. 23: end if 24: end if 25: end for Output: xms

iNEON algorithm is proposed. The algorithm is inspired by NEON [Xu et al., 2018; Allen-Zhu and Li, 2018], which is a method to find local minimum merely based on first-order oracles. Inexact NEON is a variant of NEON since the exact gradient in bilevel optimization is unavailable. However, it requires an extra nested loop to solve a subproblem that extracts a negative curvature descent direction. Besides, the gradient complexity of StocBiO with iNEON is also higher than the vanilla StocBiO. Therefore, we are motivated to propose a more efficient bilevel optimization algorithm that converges to second-order stationary point.

3 Preliminary

3.1 Notations and Assumptions

In this section we will present the notations used in this paper and introduce some basic assumptions to further illustrate the problem setting. In this paper we assume that upper-level function f(x, y) is twice differentiable. Lower-level g(x, y) is three times differentiable (only required in bilevel optimization). The partial derivative is denoted by ∇x and ∇y, e.g., ∇f(x, y) = [∇xf(x, y), ∇yf(x, y)]. Similarly, ∇2x and ∇2y represent the Hessian. ∇2xy and ∇2yx represent the Jacobian. We use ∥ · ∥2 and ∥ · ∥F to denote the spectral norm and Frobenius norm of matrix respectively. Notation ˜O(·) means the complexity after hiding logarithm terms. First, we assume that lower-level function g(x, y) is strongly-convex with respect to y so that y∗ (x) and Φ(x) can be well defined.

Assumption 1. The lower-level function g(x, y) is µstrongly-convex with respect to y, i.e., there exists a constant µ such that

for any x, y and y′ .

Notice that in minimax optimization g(x, y) is the same as -f(x, y) so we merge these two cases into one statement. With Assumption 1, objective function Φ(x) is also differentiable and the gradient is formulated as follows [Ji et al., 2021].

We can see the Hessian of g is automatically involved in the gradient of Φ. Notice that in this paper first-order method means only using the first-order information of Φ. In minimax optimization, since we always have ∇yf(x, y∗ (x)) = 0, the expression of ∇Φ(x) is simplified by

Next, we introduce the following assumptions about Lipschitz continuity of first and second order derivatives. These assumptions are commonly used in the convergence analysis of minimax and bilevel optimization [Luo et al., 2020; Luo and Chen, 2022; Ji et al., 2021; Huang et al., 2025].

Assumption 2. The gradients of component functions F(x, y; ξ) and G(x, y; ζ) are L-Lipschitz continuous, i.e., there exists a constant L such that

Assumption 3. The second order derivatives ∇2x f(x, y), ∇2xy f(x, y), ∇2y f(x, y), ∇2xy g(x, y) and ∇2y g(x, y) are ρLipschitz continuous.

The condition number κ of the hierarchical optimization problem is defined by κ = L/µ. According to previous works, in minimax optimization under Assumptions 1, 2 and 3, Φ(x) is twice differentiable. y∗ (x) is κ-Lipschitz continuous, ∇Φ(x) is LΦ-Lipschitz continuous and ∇2Φ(x) is ρΦLipschitz continuous.

According to [Ghadimi and Wang, 2018; Ji et al., 2021], we know that in bilevel optimization function y∗ (x) is also κLipschitz continuous, but we need an additional Assumptions 4 to guarantee Φ(x) has LΦ-Lipschitz gradient.

Assumption 4. The upper-level function f(x, y) is MLipschitz continuous, i.e., there exists a constant M such that

Since in this paper we study the convergence to secondorder stationary point, we also need the following Assumption 5 which is also assumed in [Huang et al., 2025] that makes function Φ(x) twice differentiable and have ρΦLipschitz Hessian. We should notice that Assumption 4 and 5 are only used for bilevel optimization.

Assumption 5. The third order derivatives ∇3xyx g, ∇3yxy g and ∇3y g are ν-Lipschitz continuous.

3.2 Relations Between Local Nash Equilibrium

In nonconvex minimax optimization, local Nash equilibrium is an important concept about the convergence criterion, which is defined as follows in [Jin et al., 2020].

Definition 1. (Local Nash equilibrium) A point (x∗ , y∗ ) is a local Nash equilibrium of Problem (2) if f(x∗ , y) ≤ f(x∗ , y∗ ) ≤ f(x, y∗ ) for any x, y that satisfy ∥x - x∗ ∥ ≤ δ and ∥y - y∗ ∥ ≤ δ.

However, as it is indicated in [Jin et al., 2020], a local Nash equilibrium may not exist in sequential games such as GAN. Hence a necessary condition of local Nash equilibrium is provided.

Definition 2. (Local minimax point) A point (x∗ , y∗ ) is a local minimax point of Problem (2) if y∗ is a local maximum of function f(x∗ , ·) and there exists a constant δ0 > 0 such that x∗ is a local minimum of function gδ(x) = max∥y-y∗ ∥≤δ f(x, y) for any 0 < δ ≤ δ0.

Next, we will prove that a saddle point of the primal function Φ(x) is not a local minimax point, which indicates the importance to find second-order stationary point of Φ(x) in minimax optimization.

Suppose x∗ is a saddle point of Φ(x) and y∗ is the maximum of f(x∗ , ·). By the definition of saddle point, we know for ∀δ > 0, there exists ∥x - x∗ ∥ < δ and Φ(x) < Φ(x∗ ). Hence for ∀δ0 > 0 and ∀δ > 0, there exists x′ such that ∥x - x∗ ∥ < min{δ0/κ, δ} and Φ(x′ ) < Φ(x∗ ). Let y′ = maxf(x′ , ·). According to Proposition 1, we have gδ0 (x′ ) = f(x′ , y′ ) = Φ(x′ ) < Φ(x∗ ) = gδ0 (x∗ ), which means x∗ is not a local minimum of function gδ0 (x). Therefore, a saddle point of the primal function Φ(x) will never be a local minimax point of Problem (2).

Proposition 1. (Lemma 4.3 in [Lin et al., 2020b]) Suppose function f satisfies Assumption 2 and Assumption 1. Then function y∗ (x) is κ-Lipschitz continuous, i.e.,

for ∀x1, x2 ∈ Rd1 . Function Φ(x) is differentiable with gradient ∇Φ(x) = ∇xf(x, y∗ (x)) and the gradient is LΦLipschitz continuous where LΦ = L + κL.

4 Algorithm for Minimax Optimization

In this section, we will propose our PRGDA algorithm for the special case of minimax optimization. The description of our PRGDA algorithm is demonstrated in Algorithm 1. Similar to SREDA, the initial value y0 is also yield by PiSARAH algorithm to make it close to y∗ (x0), which is a conventional strongly-convex optimization subproblem. In our Algorithm 2 Updater of Inner Loop (Minimax)

Input: status xt, xt-1, yt, vt-1, ut-1 and t Parameter: stepsize λ, inner loop size K, batchsize S1 and S2, period q.

• 1: Set xt,-1 = xt-1, xt,k = xt when k ≥ 0, yt,-1 = yt,0 = yt .
• 2: if mod(t, q) = 0 then
• 3: Draw S1 samples {ξ1, . . . , ξS1 }
• 4: vt,-1 = 1 S1 ∑S1 i=1 ∇xF (xt, yt; ξi),
• 5: ut,-1 = 1 S1 ∑S1 i=1 ∇yF (xt, yt; ξi).
• 6: else
• 7: vt,-1 = vt-1, ut,-1 = ut-1.
• 8: end if
• 9: for k = 0 to K -1 do
• 10: Draw S2 samples {ξ1, . . . , ξS2 }

• 12: ut,k = ut,k-1 + 1 S2 ∑S2 i=1 (∇yF(xt,k, yt,k; ξi) - ∇yF (xt,k-1, yt,k-1; ξi))
• 13: yt,k+1 = ∏Y (yt,k +λut,k).
• 14: end for
• 15: Select st = argmink ∥ ˜ Gλ(yt,k)∥.
• 16: Let yt+1 = yt,st , vt = vt,st , ut = ut,st .

Output: yt+1, vt, ut.

convergence analysis this step costs the gradient complexity of ˜O(κ2 ϵ-2 ). We use vt and ut to represent the gradient estimator of ∇xf(xt, yt) and ∇yf(xt, yt) respectively. In each iteration, yt+1, vt and ut are computed by an inner loop updater with K iterations, which is shown in Algorithm 2.

In Algorithm 2, we use the SPIDER gradient estimator to update yt,k, vt,k and ut,k. S1 is the large batchsize that is loaded every q iterations of t. S2 is the small batchsize. λ is the stepsize to update variable y. The output of the inner loop updater depends on the minimum value of the norm of ˜ Gλ(yt,k) and its corresponding index (defined in Eq. (8)). We will show that gradient estimator vt satisfies ∥vt -∇Φ(xt)∥ ≤ O(ϵ) based on this inner loop updater.

Inspired by perturbed gradient descent, our PRGDA is also composed of a descent phase and an escaping phase. In the descent phase our PRGDA algorithm follows the iterative update rule of SPIDER that xt+1 = xt - (η/∥vt∥)vt until the norm of vt satisfies ∥vt∥ ≤ O(ϵ). After the descent phase is terminated, we use ms to denote the current counter t and uniformly draw a perturbation ξ from ball B0(r) where parameter r is the perturbation radius. We add the perturbation to the current status xt and start the escaping phase. In the escaping phase, parameter tthres is maximum number of iterations of the phase and ¯D is the average moving distance which is used to determine if the escaping phase should be stopped. The stepsize of x in this phase is denoted by ηH which is typically larger than η in the descent phase. We use D to denote the accumulated squared moving distance. If the averaged squared moving distance is larger than ¯D then we pull it back (line 17 in Algorithm 1) and break the escaping phase. In this case we consider xms as a saddle point and continue to run next descent phase. Otherwise, if the escaping phase is not broken after tthres iterations, we claim that xms is a secondorder stationary point with high probability. This is because when λmin(∇2Φ(xms )) < -ϵH, the stuck region S defined by the area {ξ ∈ B0(r)| the sequence started from xms+1 = xms + ξ does not break the escaping phase} has a small volume. Informally, the stuck region S must be contained in a "narrow band" or "thin disk" in a higher dimensional sphere. Since the perturbation ξ is uniformly drawn from ball B0(r), the probability that ξ belongs to the stuck region is low.

5 Algorithm for Bilevel Optimization

In this section we propose our PRGDA algorithm to solve the more general bilevel optimization. Actually, we only need to switch the inner loop updater in Algorithm 2 to the bilevel mode, which is demonstrated in Algorithm 3. Similar to the case of minimax optimization, here we also need a initialization algorithm to initialize y0 with the cost of Gc(g, ϵ) = ˜O(κ6 ϵ-2 ) in the convergence analysis. Next we will elaborate the inner loop updater for bilevel optimization. We also use the update rule of SPIDER to compute v(1) t,k , v(2) t,k and ut,k, which represent the estimator of ∇xf(x, y), ∇yf(x, y) and ∇yg(x, y) respectively. We should notice that the large and small batchsize of computing ut,k are different from that of v(1) t,k or v(2) t,k . After the inner loop to compute yt+1, we calculate the Jacobian Jt with a batch of size S5. Then we compute vt, the estimator of ∇Φ(x) via AID. Here we follow the method used in StocBiO, which is

where Bj is the set of samples to calculate the stochastic estimator of Hessian ∇2y g(xt, yt+1).

6 Convergence Analysis

In this section we will illustrate the convergence analysis of our algorithm. First, we need to assume that Φ(x) is lower bounded by Φ∗ . Then we will present the main theorems of our PRGDA algorithm. In this paper, we set ϵH = √ρΦϵ as the tolerance of the second-order stationary point.

6.1 Main Theorem for Minimax Optimization

Theorem 1. Under Assumption 1, 2 and 3, we set stepsize η = ˜O( ϵ κL ), ηH = ˜O( 1 κL ) and λ = O( 1L ), batchsize S1 = ˜O(κ2 ϵ-2 ) and S2 = ˜O(κϵ-1 ), period q = O(ϵ-1 ), inner loop K = O(κ), perturbation radius r = min{ ˜O(√ ϵ κ3 ρ ), ˜O( ϵ κL )}, threshold tthres = ˜O( L √κρϵ ) and average movement ¯D = ˜O( ϵ2 κ2L2 ). Then our PRGDA algorithm requires ˜O(κ3 ϵ-3 ) SFO to achieve O(ϵ, √ρΦϵ) secondorder stationary point with high probability.

Algorithm 3 Updater of Inner Loop (Bilevel)

Input: status xt, xt-1, yt, v(1) t-1 , v(2) t-1 , ut-1 and t Parameter: stepsize λ and α, inner loop size K and Q, batchsize B, S1, S2, S3, S4 and S5, period q.

• 1: Set xt,-1 = xt-1, xt,k = xt when k ≥ 0, yt,-1 = yt,0 = yt .
• 2: if mod(t, q) = 0 then
• 3: Draw samples with size S1 and S3.
• 4: v(1) t,-1 = 1 S1 ∑S1 i=1 ∇xF (xt, yt; ξi),
• 5: v(2) t,-1 = 1 S1 ∑S1 i=1 ∇yF (xt, yt; ξi),
• 6: ut,-1 = 1 S3 ∑S3 i=1 ∇yG(xt, yt; ζi).
• 7: else

• 9: end if
• 10: for k = 0 to K -1 do
• 11: Draw samples with size S2 and S4.

• 14: ut,k = ut,k-1 + 1 S4 ∑S4 i=1 (∇yG(xt,k, yt,k; ζi) - ∇yG(xt,k-1, yt,k-1; ζi))
• 15: yt,k+1 = yt,k -λut,k.
• 16: end for
• 17: Select st = argmink ∥ ˜ Gλ(yt,k)∥. Let yt+1 = yt,st , v(1) t = v(1) t,st , v(2) t = v(2) t,st , ut = ut,st .
• 18: Compute Jacobian Jt = 1 ∑S5
• S5 i=1 ∇2xyG(xt, yt+1; ζi).
• 19: Compute vt via AID in Eq. (9).

Output: yt+1, vt, ut.

6.2 Main Theorem for Bilevel Optimization

Theorem 2. Under Assumption 1, 2, 3, 4 and 5, we set stepsize η = ˜O( ϵ κ3L ), ηH = ˜O( 1 κ3L ), λ = O( 1L ) and α = O( 1L ), batchsize S1 = ˜O(κ2 ϵ-2 ), S2 = ˜O(κ-1 ϵ-1 ), S3 = ˜O(κ6 ϵ-2 ), S4 = ˜O(κ3 ϵ-1 ), S5 = ˜O(κ2 ϵ-2 ) and B = ˜O(κ2 ϵ-2 ), period q = O(κ2 ϵ-1 ), inner loop K = O(κ) and Q = ˜O(κ), perturbation radius r = min{ ˜O(√ ϵ ρΦ ), ˜O( ϵ κ3L )}, threshold tthres = ˜O( κ3L √ρΦϵ ) and average movement ¯D = ˜O( ϵ2 κ6L2 ). Then our PRGDA algorithm requires complexity of Gc(f, ϵ) = ˜O(κ3 ϵ-3 ), Gc(g, ϵ) = ˜O(κ7 ϵ-3 ), JV (g, ϵ) = ˜O(κ5 ϵ-4 ) and HV (g, ϵ) = ˜O(κ6 ϵ-4 ) to achieve O(ϵ, √ρΦϵ) secondorder stationary point with high probability.

7 Experiments

In this section we conduct the matrix sensing [Bhojanapalli et al., 2016; Park et al., 2017] experiment to validate the performance of out PRGDA algorithm for solving both minimax and bilevel problem. As a result of existing study on matrix sensing problem [Ge et al., 2017], there is no spurious local minimum in this circumstance, i.e., every local minimum is a global minimum. Therefore, the capability of escaping saddle points of our algorithm can be verified by this experiment. We follow the experiment setup of [Chen et al., 2022] to recover a low-rank symmetric matrix M ∗ = U ∗ (U ∗ )T where U ∗ ∈ Rd×r . Suppose we have n sensing matrices {Ai}ni=1 with n observations bi = ⟨Ai,M ∗ ⟩. Here the inner product of two matrices is defined by the trace ⟨X, Y ⟩ = tr(XT Y ). Then the optimization problem can be defined by

7.1 Robust Optimization

Similar to the problem setting of [Yan et al., 2019], we also introduce another variable y and add a robust term to make the model robust. Therefore, the optimization problem can be formulated by

where ∆n = {y ∈ Rn |0 ≤ yi ≤ 1, ∑ni=1 yi = 1} is the simplex in Rn and Li(U) is defined in Eq. (10). Moreover, it is easy to check there is no spurious local minimum given the strict saddle property in [Ge et al., 2017].

The number of rows of matrix U is set to d = 50, d = 75 and d = 100 respectively and the number of columns is fixed as r = 3. The ground truth low-rank matrix M ∗ is generated by M ∗ = U ∗ (U ∗ )T where each entry of U ∗ is drawn from Gaussian distribution N(0, 1/d) independently. We randomly generate n = 20d samples of sensing matrices {Ai}ni=1 , Ai ∈ Rd×d from standard Gaussian distribution and calculate the corresponding labels bi = ⟨Ai,M ∗ ⟩ hence there is no noise in the synthetic data. The global minimum of loss function value Φ(U) should be 0 which can be achieved at point U = U ∗ and y = 1/n.

Following the initialization in [Chen et al., 2022], we randomly generalize a vector u0 from Gaussian distribution and multiply it by a scalar to satisfy the condition ∥u0∥ ≤ λmax(M ∗ ) where we denote λmax(·) as the maximum eigenvalue. The initial value is set to U = [u0, 0, 0]. Each optimization algorithm shares the same initialization. Apart from our PRGDA algorithm, we run three baseline algorithms, SGDA, Acc-MDA and SREDA. The code is implemented on matlab. We choose η = 0.001, ηH = 0.1, λ = 0.01, ¯D = r = 0.01, tthres = 20, K = 5, S2 = 40 and q = 25.

We evaluate the performance of each algorithm by two criteria, loss function value of Φ(U) and the ratio of distance to the optimum ∥UU T -M ∗ ∥2F /∥M ∗ ∥2F . The experimental results of these two quantities versus the number of gradient oracles are shown in Figure 1.

From the experimental results we can see SGDA, AccMDA and SREDA cannot escape saddle points because the loss function value is far away from the global minimum 0, which is equivalent to local minimum in this task because of the strict saddle property. In contrast, we can see our PRGDA algorithm eventually converges to the global optimum U ∗ and achieves the best loss function value that is close to 0, which indicates its ability to escape saddle point. Especially in the case of d = 50, we can see clearly that our PRGDA algorithm jumps out of the trap of saddle point. Besides, in our experiment we also list the smallest eigenvalue of the Hessian matrix ∇2Φ(U) for each algorithm after they have converged. The results are shown in Table 3. We can see the value λmin(∇2Φ(U)) of our method is the closest to 0 in all cases, which also verifies the performance of our PRGDA algorithm to find second-order stationary point.

7.2 Hyper-Representation Learning

We also conduct a hyper-representation learning experiment a verify the ability of our method to reach second-order stationary point in bilevel optimization. Recently, many methods in meta learning [Finn et al., 2017; Nichol et al., 2018] are designed to learn hyper-representations via two steps and separated dataset. The backbone is trained to extract better feature representations which can be applied to many different tasks. Based on these features a classifier is further learned on specific type of training data, which eventually forms a bilevel problem. In this experiment we also consider the matrix sensing task but conduct it in the hyper-representation learning manner.

The generation of U ∗ , M ∗ , Ai and bi are the same as Section 7.1. We also set d = 50, d = 75 and d = 100. The number of samples is n = 20d. We split all samples into two dataset: a train dataset D1 with 70% data and a validation dataset D2 with 30% data. We define variable x to be the first r-1 columns of U and variable y to be the last column. Then the objective function is formulated by

Here Li(·) is defined in Eq. (10) since U is the concatenation of x and y.

We follow the initialization in Section 7.1 to set x = [u0, 0] and y = 0. We compare our PRGDA algorithm with four baselines, StocBiO, MRBO, VRBO and StocBiO + iNEON. We choose η = 0.001, ηH = 0.1, λ = 0.01, ¯D = r = 0.01, tthres = 20, K = 5, S2 = 40 and q = 25. We also use the ration of distance to optimum, i.e. ∥UU T -M ∗ ∥2F /∥M ∗ ∥2F as the metric to evaluate the performance. The experimental results are shown in Figure 2.

From the experimental results we can see our PRGDA algorithm shows the best performance to reach second-order stationary point and approach the expected optimum. MRBO and VRBO do not escape saddle points during the experiment. In the case of d = 50, StocBiO performs better than MRBO and VRBO because the randomness of stochastic gradient serves as a kind of perturbation, while in variancereduced algorithms the gradient estimator is closer to the full gradient. This result indicates the necessity of our method to make variance-reduced bilevel algorithm escape saddle points. StocBiO + iNEON also escapes saddle point probably but its convergence is slower than our method.

8 Conclusion

In this paper, we investigate the methods to achieve secondorder optimality for nonconvex minimax and bilevel optimization. We propose a new algorithm PRGDA for stochastic nonconvex hierarchical optimization which is the first algorithm to find second-order stationary point for stochastic nonconvex-strongly-concave minimax optimization. In nonconvex-strongly-convex bilevel optimization, our method also achieves better gradient complexity to find local minimum. We prove that our method obtains the gradient complexity of ˜O(ϵ-3 ) to achieve O(ϵ, √ρΦϵ) second-order stationary point, which matches the best results of searching first-order stationary point under same conditions. We also conduct two numerical experiments to verify the performance of our algorithm.

Acknowledgements

This work was partially supported by NSF IIS 2347592, 2348169, DBI 2405416, CCF 2348306, CNS 2347617.