The study of saddle flow dynamics traces back to the seminal work [1] and contributes to the fundamentals of mathematical optimization from a dynamical system perspective. Saddle flow dynamics aim to search for minimax saddle points by combining gradient descent with gradient ascent on two respective subsets of variables (that form a partition). This approach is the basis of primal-dual methods used for solving constrained optimization problems [2] and bestresponse dynamics used for locating Nash equilibria in zerosum games [3], [4]. As a result, saddle flow dynamics have been widely used for resources allocation and stabilizing controller design in a variety of areas, including power systems [5], [6], communication networks [7], [8], and cloud computing [9].
To fully understand the dynamic behavior of saddle flows (and its variants), the convergence properties have been extensively studied. Preliminary results [10]- [13] are centered on asymptotic stability, primarily using advanced analytical tools to provide insights into the important special case of primal-dual dynamics. More recent efforts [14]- [16] have been made to study the exponential convergence of saddle flow dynamics for not only theory development but also practical uses. For example, [14], [15] have shown particular exponential convergence in the absence of strong convexity. [2] has further extended the study to discrete-time problems, showing linear convergence -the counterpart of exponential convergence in continuous time.
A fundamental challenge in studying these exponentially convergent dynamics is to estimate their decay rates using several algorithmic variants [17]- [23]. A vast amount of literature resorts to proximal regularization, especially when handling non-smooth functions. In particular, [17] only shows the existence of exponential convergence rates related to strong convexity but fails to provide an estimate. Lower bounds on the convergence rate have been developed in [18], [19] and [20], respectively, using saddle flow variants and frequency-domain Integral Quadratic Constraint (IQC) approaches. [21] further establishes a best-case upper-bound rate with a Lyapunov method in the time domain. Another set of results [22], [23] develop a novel projection on the standard Lagrangian and derive lower bounds on the decay rates when proving (semi-)global exponential stability for the augmented primal-dual dynamics. By and large, the rate estimates provided in the literature depend on strong convexity, regularization coefficients, and singular values of any coupling matrices. Notably, the constant of strong convexity seems to be a bottleneck inherent in the convergence rate of saddle flow dynamics.
In this paper, we focus particularly on a class of bilinearly coupled minimax problems with strongly convexlinear objective functions and design an algorithm based on carefully designed saddle flow dynamics that exploits the problem structure and achieves an exponential convergence rate beyond the strong convexity constant. The algorithm involves two sequential steps -change-of-variable conditioning and proximal regularization -that jointly enhance the strong convexity-strong concavity of a given objective function. Building upon our recent results [24], this property immediately results in a lower-bound estimate for the convergence rate that breaks the stereotype limit. In addition, such a performance guarantee makes the algorithm suited for a class of weakly convex-weakly concave functions with exponential convergence to a stationary point. Extensive simulations are run to compare our algorithm with existing methods and the actual rate with its lower-bound rate estimate.
Our contributions are summarized as follows. (i) For bilinearly coupled minimax problems with strongly convex-linear objective functions, we design a saddle flow-based algorithm with an exponential convergence rate beyond the given strong convexity constant. (ii) We provide an explicit lower-bound estimate for the exponential convergence rate of the algorithm, which offers a guideline for choosing parameters to further optimize the rate. (iii) The accelerated rate guarantee accommodates weakly convex-weakly concave functions, thus showing exponential convergence in nonconvex-nonconcave scenarios. (iv) Numerical results validate the proposed algorithm's superiority to many existing methods and further show that the lower-bound rate estimate is almost tight. Notation: Let R be the set of real numbers. I n ∈ R n×n denotes the identity matrix of size n. Given a twice differentiable function L(x, y) ∈ C 2 with L : R n × R m → R, we use ∂ ∂x L(x, y) ∈ R 1×n and ∂ ∂y L(x, y) ∈ R 1×m to denote the partial derivatives with respect to x and y, respectively. We further define ∇ x L(x, y) := ∂ ∂x L(x, y) T . Meanwhile, ∂ 2
∂x 2 L(x, y) ∈ R n×n and ∂ 2 ∂y 2 L(x, y) ∈ R m×m represent the second-order partial derivatives of L(x, y) with respect to x and y, respectively.
In this paper, we aim to solve bilinearly coupled minimax problems of the form:
where η > 0 denotes a constant, b ∈ R m is a constant vector, and A ∈ R m×n is the coupling matrix. f (x) : R n → R is assumed to be twice differentiable. We assume a saddle point (x ⋆ , y ⋆ ) exists. The minimax problem (1) is equivalent to an affine equality-constrained optimization problem:
where y ∈ R m represents the dual variable. Note that L(x, y) is the (weighted) Lagrangian of the problem (2) and (x ⋆ , ηy ⋆ ) is an optimal primal-dual solution. The primary goal is to locate a saddle (min-max) point of L(x, y), as indicated in (1). We will particularly focus on saddle flow dynamics that run a continuous-time version of gradient descent ascent on L(x, y): Definition 1 (Saddle Flow Dynamics). The saddle flow dynamics on L(x, y) are defined as the following dynamic law:
In the following subsections, we will develop an algorithm based on saddle flow dynamics that exponentially converges to a saddle point of L(x, y).
Our algorithm is defined on a set of new auxiliary variables. Three auxiliary variables u ∈ R n , p, v ∈ R m are introduced, and their connection with the original variables x and y is briefly illustrated in Figure 1.
, . First, we propose a change of variables to transform (x, y) into (u, p) dictated by:
where α > 0 is a positive constant, and B ∈ R m×m is a given full-rank transformation matrix, with σB I ⪰ BB T ⪰ σ B I for some σB > 0 and σ B > 0. We obtain L(u, p) as the counterpart of the Lagrangian L(x, y) in the (u, p) space:
This conditioning step exploits the linearity in y. Second, we introduce a proximal regularizer on the variable p to obtain a Moreau envelope function defined on (u, v) in the following form:
where ρ > 0 denotes a constant regularization coefficient. The saddle flow dynamics on the Moreau envelope are expected to reach an equilibrium point (u ⋆ , v ⋆ ), which immediately implies a saddle point (x ⋆ , y ⋆ ) of the original Lagrangian L(x, y) via the inverse change of variables.
More specifically, we define our algorithm as the following dynamic law on (u, v):
where p is a shorthand for the mapping p(u, v) : R n ×R m → R m that represents the unique maximizer in (6) (uniqueness to be shown later) subject to
We will provide the conditions under which the algorithm ( 7) is guaranteed to achieve the goal and analyze its performance in the next subsection.
Before stating the main result, we outline the assumptions required for the algorithm (7) that basically characterize the objective function f (x) and the matrix A of the affine equality constraint in (2).
Given Assumption 1, L(x, y) is a strongly convex-linear function.
Assumption 2. The coupling matrix A is full row rank with σA I ⪰ AA T ⪰ σ A I for some σA > 0 and σ A > 0.
Given Assumptions 1 and 2, the saddle point (x ⋆ , y ⋆ ) of L(x, y) is unique since the problem (2) meets the linear independence constraint qualification and admits a unique primal-dual solution [25], [26]. Note that the saddle point (x ⋆ , y ⋆ ) of L(x, y) is also the equilibrium point of the corresponding saddle flow dynamics (3).
Our first argument is that an equilibrium point (u ⋆ , v ⋆ ) of the algorithm ( 7) is immediately related to the target saddle point (x ⋆ , y ⋆ ), as the following theorem suggests.
Theorem 1 (Characterization of Saddle Point). Let Assumptions 1 and 2 hold. Given an equilibrium point (u ⋆ , v ⋆ ) of the dynamic law (7), the unique saddle point (x ⋆ , y ⋆ ) of L(x, y) can be attained by
The proof of Theorem 1 is provided in Appendix A, which is an immediate consequence of the algorithm design in Section III.
Given Theorem 1, it remains to show the convergence of the algorithm (7) to an equilibrium point, which we establish along with a lower-bound rate estimate in the following theorem.
Theorem 2 (Convergence). Let Assumptions 1 and 2 hold. Given any positive constants η and α that satisfy 2η > αl f , the dynamic law (7) globally exponentially converges to a unique equilibrium point (u ⋆ , v ⋆ ). More precisely, given w := (u, v),
holds with the lower-bound rate
where γ 1 is a constant with
The proof will become clear as we develop the algorithm (7) in detail in Section III.
Theorem 2 not only guarantees the exponential convergence of the algorithm (7), but also implies the following corollary that suggests an accelerated rate beyond the constant of the strong convexity of f (x).
Corollary 3. Given any positive constants ρ, η, σB and σ B that satisfy (7) achieves a convergence rate lower bounded by c > m f .
In Corollary 3, the first condition ensures that the righthand side of the inequality in the second condition is positive. Technically, the second condition can be derived from the expression for the rate c. The intuition will be provided later in Section III.
Remark 1. The parameters that satisfy all the conditions always exist. A straightforward way is to pick a sufficiently small σB and a sufficiently large ρ. For example, the following set of parameters
satisfy all the conditions.
In this section, we explain the development of the algorithm, which naturally provides insights into the accelerated convergence performance of the algorithm (7). Inspired by [27] and [28], we aim to enhance the convexity-concavity of the original Lagrangian in (1) by exploiting the bilinear couplings.
Recall the change of variables (4) in the first step. We obtain the expression for the Lagrangian L(x, y) in the (u, p) space as L(u, p) in (5). The role of the change-of-variable conditioning (4) is captured below. Lemma 4. Let Assumptions 1 and 2 hold. Given any positive constants η and α that satisfy 2η > αl f , L(u, p) is m fstrongly convex in u and
-strongly concave in p.
Proof. Lemma 4 follows immediately from the calculation of the respective second-order partial derivatives:
■ Lemma 4 suggests that the change-of-variable conditioning enhances the concavity (from linearity) while maintaining the strong convexity for the objective function. In fact, the constant of strong concavity of L(u, p) in p could be further optimized with proper choices of parameters. See an example below.
Remark 2. When first fixing η and σB , setting α = η/l f gives the largest possible constant of strong concavity
Furthermore, L(u, p) is relevant since its unique saddle point (u ⋆ , p ⋆ ) is directly connected with the saddle point (x ⋆ , y ⋆ ) of the original Lagrangian L(x, y), as the corollary states.
Corollary 5. (u ⋆ , p ⋆ ) is the unique saddle point of L(u, p) if and only if (x ⋆ , y ⋆ ) is the unique saddle point of L(x, y) with
This subsection builds on L(u, p) and aims to further improve its strong convexity-strong concavity. As the second step, we apply proximal regularization on p to obtain the Moreau envelope of L(u, p):
where p is the solution to the optimality condition (8). Note that p is unique due to the strong concavity of L(u, p) in p.
It can be efficiently computed, either using its closed form when it is available, i.e., the inverse function of ∇f (•) has an explicit expression, or using numerical methods since ∇f (•) is strictly monotone. The following lemma characterizes the impact of proximal regularization on the convexity-concavity property. Lemma 6. Let Assumptions 1 and 2 hold. Given any positive constants η and α that satisfy 2η > αl f , L(u, v) is c 1strongly convex in u and c 2 -strongly concave in v with
Proof. Since L(x, y) is twice differentiable, we are allowed to derive the partial derivatives of L(u, v) as
Using Danskin's Theorem, we can further compute the second-order partial derivatives. Define
which is a symmetric matrix. Then we obtain
where J u p and J v p are the partial derivatives of p(u, v) with respect to u and v, respectively. By taking partial derivatives of both sides of the optimality condition (8), we have
The explicit second-order partial derivatives from above allow us to characterize the strong convexity-strong concavity of L(u, v):
Compared with the strong convexity-strong concavity of L(u, p) identified in Lemma 4, Lemma 6 reveals an important role of proximal regularization here: it increases the constant of strong convexity while decreasing the constant of strong concavity.
Similarly, we can establish the immediate connection between the saddle points of L(u, p) and L(u, v) as follows.
Note that given Lemmas 4 and 6 as well as Corollaries 5 and 7, Theorem 1 becomes straightforward. We summarize the brief proof in Appendix A.
Now it becomes clear that the algorithm ( 7) is essentially the saddle flow dynamics on L(u, v). The exponential convergence of the algorithm established in Theorem 2 could be built upon our recent result [24], which we briefly present here for completeness.
Lemma 8. Given a twice differentiable function S(x, y) with m-strong convexity in x and q-strong concavity in y, the saddle flow dynamics on S(x, y):
with z := (x, y), are globally exponentially stable. More precisely, ∥z(t)z ⋆ ∥ ≤ ∥z(0)z ⋆ ∥e -rt holds with the lower-bound rate r := min{m, q}.
The proof of Lemma 8 is available in [24].
In virtue of Lemma 8, the strong convexity-strong concavity of L(u, v) in Lemma 6 immediately implies Theorem 2.
Rate Optimization: According to the parameter conditions in Corollary 3, there exist two constants 0
Then the strong convexity and strong concavity constants of L(u, v) can be rewritten as
where κ A := σA /σ A is the condition number of AA T . Since c 1 decreases while c 2 increases, both monotonically in ρ, and there is bound to be one crossing point (given the respective ranges), we can obtain the fastest convergence rate in terms of ρ by setting c 1 = c 2 . It will be achieved at ρ ⋆ that satisfies
with L =: 2ηαα 2 l and M := 2ηαα 2 m f . Further optimization over the other parameters is intricate and will be left for future studies.
In this section, we study a broader class of bilinearly coupled minimax problems:
where A ∈ R m×n is the coupling matrix while f (x) : R n → R and g(y) : R m → R are assumed to be twice differentiable, but not necessarily convex. We assume a stationary point exists.
Searching for saddle/stationary points of nonconvexnonconcave functions is a popular goal in the machine learning community where such objective functions are common [29], [30]. While in general locating or even approximating stationary points is hard, certain special structures of objective functions, e.g., bilinear couplings, make this problem still tractable [28], [31], [32].
In this section, we extend the results of this paper to solve the problem (19), in the case where L(x, y) is weakly convex-weakly concave and smooth. Assumption 3. f (x) ∈ C 2 is m f -weakly convex and l fsmooth, i.e., l f I ⪰ ∇ 2 f (x) ⪰ -m f I with l f > 0 and m f > 0. Similarly, g(y) ∈ C 2 is m g -weakly convex and l g -smooth, i.e., l g I ⪰ ∇ 2 g(y) ⪰ -m g I with l g > 0 and m g > 0.
Remark 3. Note that, the objective function L(x, y) is convex-concave when m f = m g = 0 holds and strongly convex-strongly concave when m f < 0 and m g < 0 hold. This study specifically addresses the scenarios in which m f > 0 and m g > 0 hold, thereby quantifying the degree of nonconvexity-nonconcavity of L(x, y).
In what follows we use a pipeline similar to that in Sections II and III, highlighting mostly the main changes due to Assumption 3. As before, applying the change of variables (4), we obtain the corresponding objective function L(u, p) in the (u, p) space:
Even though in general L(u, p) is not necessarily concave in p, one can compute the second-order partial derivatives (similar to the procedure to prove Lemma 4) and show that L(u, p) is
-strongly concave in p whenever σ A (2α-α 2 l f )-m g > 0 holds. This guarantees the existence and uniqueness of the following Moreau envelop:
where p⋆ , a shorthand for p⋆ (u, v), is the unique maximizer in (21). We will show shortly that under mild conditions, L(u, v) in ( 21) is strongly convex-strongly concave, and thus the corresponding saddle flow dynamics are guaranteed to convergence to a unique stationary/saddle point (u ⋆ , v ⋆ ). The following theorem relates such a point with a stationary point of L(x, y).
is a stationary point of L(u, p). Moreover, (x ⋆ , y ⋆ ), as described by (9), is a stationary point of L(x, y).
The proof is provided in Appendix B. Theorem 9 guarantees that if the saddle flow dynamics on L(u, v), i.e.,
converge, then the corresponding (x ⋆ , y ⋆ ) will be a stationary point of L(x, y).
0.0 0.5 1.0 1.5 2.0 2.5 3.0 Time -150 -100 -50 0 50 Variables x 1 x 2 x 3 x 4 x 5 0.50 0.75 1.00 1.25 1.50 -0.2 0.0 0.2 (a) Trajectories of part of variables. 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Time 10 -8 10 -6 10 -4
10 -2 10 0 Distance to Equilirium Point Proximal Method PDGD Strong Convexity mf Our Estimate Rate Our Method (b) Convergence rate.
Fig. 2: Results for a strongly convex-linear minimax problem.
We finalize this section by providing a convergence guarantee for (22).
Theorem 10. Let Assumptions 2 and 3 hold. Given any positive constants α, ρ, σB and σ B that satisfy 1)
then the dynamic law (22) exponentially converges to a unique equilibrium point (u ⋆ , v ⋆ ). More precisely, given w := (u, v), ∥w(t)w ⋆ ∥ ≤ ∥w(0)w ⋆ ∥e -ct holds with the lower-bound rate
The proof is similar to that of Theorem 2.
Remark 4. Since the original objective function L(x, y) is nonconvex-nonconcave, there might be multiple stationary points. However, L(u, v) becomes strongly convex-strongly concave and has a unique saddle point. As a result, our algorithm only obtains one particular stationary point of L(x, y).
In this section, we numerically validate our algorithm (7) to solve the strongly convex-linear minimax problem (1) and compare it with some existing algorithms in the literature. Precisely, we set n = 5 and m = 4, and adopt a quadratic objective function f (x) = 1 2 x T Qx. Here Q ∈ R 5×5 is given by Q = 5I 5 + Q T 0 Q 0 , and Q 0 ∈ R 5×5 is generated from a Gaussian random matrix. Similarly, A ∈ R 4×5 and b ∈ R 4 are also randomly generated. In this particular test run, we have the strong convexity constant m f = 0.5 and the smoothness constant l f = 0.693. The rest of the parameters are set according to (10) and the initial point is chosen arbitrarily.
The convergence results are shown in Fig. 2. It can be observed that the trajectories of the saddle flow dynamics converge rapidly to an equilibrium point in Fig. 2(a). The distance to the equilibrium point exponentially decays, as shown in Fig. 2(b). Our lower-bound estimate is close to the actual rate and almost tight. Moreover, our method is superior to many existing algorithms, including the proximal method [20] and the primal-dual gradient dynamics (PDGD) in [23], the convergence rates of which are limited by the strong convexity constant m f .
This paper studies bilinearly coupled minimax problems and proposes an accelerated algorithm to locate saddle points for strongly convex-linear objective functions. The algorithm implements saddle flow dynamics on a transformed objective function, obtained from two sequential steps: change-ofvariable conditioning and proximal regularization. The first step exploits the bilinear couplings to enhance the concavity while maintaining the strong convexity for the objective function. The second step then shifts some of the strong concavity to enhance the strong convexity. We propose the conditions under which our algorithm can achieve a convergence rate beyond the stereotype limit (the strong convexity constant). The same design strategy can also be applied to construct a stationary point searching algorithm for a general class of weakly convex-weakly concave functions.